Voyager II¶
In [2]:
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
import math as m
from scipy.optimize import fsolve, brentq, newton, least_squares, fmin
from matplotlib.lines import Line2D
from matplotlib import colors
from matplotlib import cm
from matplotlib.colors import Normalize
from scipy.interpolate import griddata
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.ticker import ScalarFormatter
import time
import warnings
from tqdm import tqdm
import pickle
import gc
from astropy import units as u
from astropy.time import Time
from astropy.coordinates import solar_system_ephemeris
solar_system_ephemeris.set("jpl")
from poliastro.bodies import *
from poliastro.ephem import Ephem
from poliastro.frames import Planes
from poliastro.maneuver import Maneuver
from poliastro.twobody import Orbit
from poliastro.util import norm, time_range
from poliastro.plotting.static import StaticOrbitPlotter
from poliastro.core.elements import rv2coe, coe2rv
from datetime import datetime
from itertools import permutations, combinations
from random import shuffle
import random
import statistics
import seaborn as sns
from IPython.display import HTML
import matplotlib.animation as animation
from poliastro.core.iod import vallado
from poliastro.core.elements import rv2coe
import datetime
import os
from pathlib import Path
from IPython.display import HTML, display
from matplotlib.animation import PillowWriter
import io
from base64 import b64encode
def show(variable, units):
def get_var_name(var):
for name, value in globals().items():
if value is var:
return name
print(get_var_name(variable), "=", variable, units)
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
In [3]:
def set_font_sizes(legend_size=14):
# Set LaTeX for text rendering
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
# Set font sizes based on legend size
plt.rc("legend", fontsize=legend_size) # legend fontsize (base size)
plt.rc(
"axes", titlesize=legend_size * 1.2
) # fontsize of the axes title (20% larger)
plt.rc(
"axes", labelsize=legend_size * 1.0
) # fontsize of the x and y labels (same as legend)
plt.rc(
"xtick", labelsize=legend_size * 0.8
) # fontsize of the tick labels (20% smaller)
plt.rc(
"ytick", labelsize=legend_size * 0.8
) # fontsize of the tick labels (20% smaller)
plt.rc(
"figure", titlesize=legend_size * 1.4
) # fontsize of the figure title (40% larger)
annotation_size = (
legend_size * 0.9
) # annotation text size (10% smaller than legend)
plt.rc("font", size=annotation_size) # default font size for annotations
# Return the settings dict for reference (optional)
settings = {
"legend": legend_size,
"axes.title": legend_size * 1.2,
"axes.label": legend_size * 1.0,
"xtick.label": legend_size * 0.8,
"ytick.label": legend_size * 0.8,
"figure.title": legend_size * 1.4,
"annotation": annotation_size,
}
return settings
In [4]:
# Create an output directory handler
class OutputManager:
def __init__(self, base_dir="V1"):
# Initialise the output manager with a base directory
self.base_dir = Path(base_dir)
self.base_dir.mkdir(exist_ok=True)
print(f"Output directory set to: {self.base_dir.absolute()}")
def save_file(self, filename, dpi=500, bbox_inches="tight", **kwargs):
# Save a figure with consistent settings.
# Ensure the filename has an extension
if not any(filename.endswith(ext) for ext in [".png", ".jpg", ".pdf", ".svg"]):
filename = f"{filename}.png"
# Create the full path
full_path = self.base_dir / filename
# Save the figure
plt.savefig(full_path, dpi=dpi, bbox_inches=bbox_inches, **kwargs)
print(f"Figure saved to: {full_path}")
return full_path
def get_path(self, filename):
"""Get the full path for a file."""
return self.base_dir / filename
# Create an instance of the output manager
output_figures = OutputManager("Voyager 2 - Figures")
output_pickle = OutputManager("Raw Results")
output_ani = OutputManager("Voyager 2 - Animations")
Output directory set to: /Users/snath2/Desktop/Aero & Astro/Year 3/IP/Final IP/Optimisation Algorithms/Voyager 2 - Figures Output directory set to: /Users/snath2/Desktop/Aero & Astro/Year 3/IP/Final IP/Optimisation Algorithms/Raw Results Output directory set to: /Users/snath2/Desktop/Aero & Astro/Year 3/IP/Final IP/Optimisation Algorithms/Voyager 2 - Animations
In [5]:
global_flag = False
In [6]:
testing_flag = False
In [7]:
mu_earth = Earth.k.to(u.km**3 / u.s**2).value
mu_jupiter = Jupiter.k.to(u.km**3 / u.s**2).value
mu_saturn = Saturn.k.to(u.km**3 / u.s**2).value
mu_uranus = Uranus.k.to(u.km**3 / u.s**2).value
mu_neptune = Neptune.k.to(u.km**3 / u.s**2).value
mu_sun = Sun.k.to(u.km**3 / u.s**2).value
mu_sun
Out[7]:
132712442099.00002
In [8]:
mu = mu_sun # Sun's gravitational parameter (km³/s²)
mu
Out[8]:
132712442099.00002
In [9]:
# Define parula colormap (MATLAB-style) with an option to reverse
def parula_map(reverse=False):
# RGB values
cm_data = [
[0.2081, 0.1663, 0.5292],
[0.2116, 0.1898, 0.5777],
[0.2123, 0.2138, 0.6270],
[0.2081, 0.2386, 0.6771],
[0.1959, 0.2645, 0.7279],
[0.1707, 0.2919, 0.7792],
[0.1253, 0.3242, 0.8303],
[0.0591, 0.3598, 0.8683],
[0.0117, 0.3875, 0.8820],
[0.0060, 0.4086, 0.8828],
[0.0165, 0.4266, 0.8786],
[0.0329, 0.4430, 0.8720],
[0.0498, 0.4586, 0.8641],
[0.0629, 0.4737, 0.8554],
[0.0723, 0.4887, 0.8467],
[0.0779, 0.5040, 0.8384],
[0.0793, 0.5200, 0.8312],
[0.0749, 0.5375, 0.8263],
[0.0641, 0.5570, 0.8240],
[0.0488, 0.5772, 0.8228],
[0.0343, 0.5966, 0.8199],
[0.0265, 0.6137, 0.8135],
[0.0239, 0.6287, 0.8038],
[0.0231, 0.6418, 0.7913],
[0.0228, 0.6535, 0.7768],
[0.0267, 0.6642, 0.7607],
[0.0384, 0.6743, 0.7436],
[0.0590, 0.6838, 0.7254],
[0.0843, 0.6928, 0.7062],
[0.1133, 0.7015, 0.6859],
[0.1453, 0.7098, 0.6646],
[0.1801, 0.7177, 0.6424],
[0.2178, 0.7250, 0.6193],
[0.2586, 0.7317, 0.5954],
[0.3022, 0.7376, 0.5712],
[0.3482, 0.7424, 0.5473],
[0.3953, 0.7459, 0.5244],
[0.4420, 0.7481, 0.5033],
[0.4871, 0.7491, 0.4840],
[0.5300, 0.7491, 0.4661],
[0.5709, 0.7485, 0.4494],
[0.6099, 0.7473, 0.4337],
[0.6473, 0.7456, 0.4188],
[0.6834, 0.7435, 0.4044],
[0.7184, 0.7411, 0.3905],
[0.7525, 0.7384, 0.3768],
[0.7858, 0.7356, 0.3633],
[0.8185, 0.7327, 0.3498],
[0.8507, 0.7299, 0.3360],
[0.8824, 0.7274, 0.3217],
[0.9139, 0.7258, 0.3063],
[0.9450, 0.7261, 0.2886],
[0.9739, 0.7314, 0.2666],
[0.9938, 0.7455, 0.2403],
[0.9990, 0.7653, 0.2164],
[0.9955, 0.7861, 0.1967],
[0.9880, 0.8066, 0.1794],
[0.9789, 0.8271, 0.1633],
[0.9697, 0.8481, 0.1475],
[0.9626, 0.8705, 0.1309],
[0.9589, 0.8949, 0.1132],
[0.9598, 0.9218, 0.0948],
[0.9661, 0.9514, 0.0755],
[0.9763, 0.9831, 0.0538],
]
if reverse:
cm_data = cm_data[::-1] # Reverse the list if specified
return colors.ListedColormap(cm_data)
1) Choose the epochs $T_1$ to $T_N$, corresponding to all the planets $P_1$ to $P_N$ in the mission. This is decided later in the optimisation algorithm.¶
In [11]:
# Voyager II real mission timing
scale = "tt"
date_departure = "1977-08-23 11:29:11"
date_flyby_jupiter = "1979-07-09 22:29:51"
date_flyby_saturn = "1981-08-26 03:24:57"
date_flyby_uranus = "1986-08-24 17:59:47"
date_arrival = date_flyby_neptune = "1989-08-25 03:56:36"
1.1) Convert epochs to MJD¶
In [13]:
from astropy.time import Time
# e.g T_et = "1979-03-05 12:05:26"
def convert_ET(T_et):
# Defines time scale for astropy.time of ET, which approx equal to TT
T_tt = Time(T_et, scale="tt")
# Calculate TDB (only 1-2 milliseconds different from ET/TT)
T_tdb = T_tt.tdb
# Calculates the TAI (T_tai = T_tt - 32.184s)
T_tai = T_tt.tai
# Calculates the JD from epoch 1st January 4713 BC at 12:00 TT
T_jd = T_tt.jd
# Calculates MJD
T_mjd = T_jd - 2400000.5
# Calculate UTC using NASA's Sprice Time Subsystem
T_utc = T_tt.utc
return {"TDB": T_tdb, "TAI": T_tai, "JD": T_jd, "MJD": T_mjd, "UTC": T_utc}
scale = "mjd"
date_departure_mjd = convert_ET(date_departure)["MJD"]
date_flyby_jupiter_mjd = convert_ET(date_flyby_jupiter)["MJD"]
date_flyby_saturn_mjd = convert_ET(date_flyby_saturn)["MJD"]
date_flyby_uranus_mjd = convert_ET(date_flyby_uranus)["MJD"]
date_flyby_neptune_mjd = convert_ET(date_flyby_neptune)["MJD"]
show(date_departure_mjd, "MJD")
show(date_flyby_jupiter_mjd, "MJD")
show(date_flyby_saturn_mjd, "MJD")
show(date_flyby_uranus_mjd, "MJD")
show(date_flyby_neptune_mjd, "MJD")
date_departure_mjd = 43378.478599537164 MJD date_flyby_jupiter_mjd = 44063.93739583343 MJD date_flyby_saturn_mjd = 44842.142326388974 MJD date_flyby_uranus_mjd = 46666.74984953692 MJD date_flyby_neptune_mjd = 47763.164305555634 MJD
1.2) Find the $\Delta t$ for each mission trajectory¶
In [15]:
def find_tof(start_date_days, end_date_days):
delta_t = (end_date_days - start_date_days) * (60 * 60 * 24)
return delta_t
delta_t_earth_jupiter = find_tof(date_departure_mjd, date_flyby_jupiter_mjd)
delta_t_jupiter_saturn = find_tof(date_flyby_jupiter_mjd, date_flyby_saturn_mjd)
delta_t_saturn_uranus = find_tof(date_flyby_saturn_mjd, date_flyby_uranus_mjd)
delta_t_uranus_neptune = find_tof(date_flyby_uranus_mjd, date_flyby_neptune_mjd)
print("Earth-Jupiter:", delta_t_earth_jupiter, "s")
print("Jupiter-Saturn:", delta_t_jupiter_saturn, "s")
print("Saturn-Uranus:", delta_t_saturn_uranus, "s")
print("Uranus-Neptune:", delta_t_uranus_neptune, "s")
Earth-Jupiter: 59223639.999997616 s Jupiter-Saturn: 67236905.99999875 s Saturn-Uranus: 157646089.9999827 s Uranus-Neptune: 94730209.00001675 s
2) Using the known epochs $T_1$ to $T_N$, calculate their respective planetary state vectors $\{\vec{\mathbf{r}},\vec{\mathbf{V}}^{\,(\text{P})}\}$ using Algorithm 3.3.¶
In [17]:
def get_planet_state_vector(planet, date_time, scale):
# Convert the date_time in MJD to modern TDB time scale
date = Time(date_time, format="mjd", scale=scale).tdb
# Generate the planet's ephemeris at the given time
planet_ephem = Ephem.from_body(planet, date, plane=Planes.EARTH_ECLIPTIC)
# Extract the position and velocity vectors
r, v = planet_ephem.rv(date)
r = r.to(u.km)
v = v.to(u.km / u.s)
return r, v
scale = "tt"
r_earth, v_earth = get_planet_state_vector(Earth, date_departure_mjd, scale)
r_jupiter, v_jupiter = get_planet_state_vector(Jupiter, date_flyby_jupiter_mjd, scale)
r_saturn, v_saturn = get_planet_state_vector(Saturn, date_flyby_saturn_mjd, scale)
r_uranus, v_uranus = get_planet_state_vector(Uranus, date_flyby_uranus_mjd, scale)
r_neptune, v_neptune = get_planet_state_vector(Neptune, date_flyby_neptune_mjd, scale)
In [18]:
print("Earth:")
display(r_earth, v_earth)
print("Jupiter:")
display(r_jupiter, v_jupiter)
print("Saturn:")
display(r_saturn, v_saturn)
print("Uranus:")
display(r_uranus, v_uranus)
print("Neptune:")
display(r_neptune, v_neptune)
Earth:
$[1.3203063 \times 10^{8},~-75056544,~-12567.941] \; \mathrm{km}$
$[14.170108,~25.826829,~0.00044418834] \; \mathrm{\frac{km}{s}}$
Jupiter:
$[-5.8784089 \times 10^{8},~5.3724892 \times 10^{8},~10951986] \; \mathrm{km}$
$[-8.9623414,~-9.0425799,~0.23792779] \; \mathrm{\frac{km}{s}}$
Saturn:
$[-1.4038366 \times 10^{9},~-2.9149948 \times 10^{8},~60883020] \; \mathrm{km}$
$[1.4494725,~-9.4742958,~0.10784132] \; \mathrm{\frac{km}{s}}$
Uranus:
$[-4.2717244 \times 10^{8},~-2.8335219 \times 10^{9},~-4983585.4] \; \mathrm{km}$
$[6.6827976,~-1.3326801,~-0.091621069] \; \mathrm{\frac{km}{s}}$
Neptune:
$[8.9742641 \times 10^{8},~-4.4292466 \times 10^{9},~70524745] \; \mathrm{km}$
$[5.2919942,~1.1093555,~-0.14477766] \; \mathrm{\frac{km}{s}}$
In [19]:
r_earth, v_earth = r_earth.value, v_earth.value
r_jupiter, v_jupiter = r_jupiter.value, v_jupiter.value
r_saturn, v_saturn = r_saturn.value, v_saturn.value
r_uranus, v_uranus = r_uranus.value, v_uranus.value
r_neptune, v_neptune = r_neptune.value, v_neptune.value
3) Now, to define the transfer trajectory we separate them into legs: $P_1\text{-}P_2$, $P_2\text{-}P_3$, ..., $P_{N-1}\text{-}P_N$. For each leg, we can use the solution to Lambert's Problem (Algorithm 3.4) to get its transfer velocity vectors. This fully defines each leg, where for instance $P_1\text{-}P_2$ leg, we have the following parameters $\big\{\vec{\mathbf{r}}_{1}, \vec{\mathbf{r}}_{2}, \vec{\mathbf{V}}^{\,(\text{P})}_1,\vec{\mathbf{V}}^{\,(\text{P})}_2, (\vec{\mathbf{v}})_{D}, (\vec{\mathbf{v}}_{T_2})_{A}, \Delta t\big\}$.¶
https://docs.poliastro.space/en/stable/autoapi/poliastro/core/iod/index.html
In [21]:
# from poliastro.core.iod import vallado
# # Lambert Solver
# def lambert(r1, r2, delta_t):
# v_d, v_a = vallado(mu_sun, r1, r2, delta_t, M=0, prograde=True, lowpath=True, numiter=1e6, rtol=1e-8)
# return v_d, v_a
import numpy as np
from scipy.optimize import brentq, newton, least_squares
# Known parameters [r1, r2, delta_t]
mu = 1.327e11 # Sun's gravitational parameter (km^3/s^2)
# Stumpff functions
def C(z):
if z > 0:
return (1 - np.cos(np.sqrt(z))) / z
elif z < 0:
return (np.cosh(np.sqrt(-z)) - 1) / (-z)
else:
return 0.5
def S(z):
if z > 0:
return (np.sqrt(z) - np.sin(np.sqrt(z))) / (z**1.5)
elif z < 0:
return (np.sinh(np.sqrt(-z)) - np.sqrt(-z)) / ((-z) ** 1.5)
else:
return 1 / 6
# Compute y parameter
def compute_y(z, r1, r2, A):
C_z = max(C(z), 1e-8) # Ensure no division by zero
return max(r1 + r2 + A * (z * S(z) - 1) / np.sqrt(C_z), 1e-8)
# Time-of-flight equation
def F_lambert(z, delta_t, r1, r2, A):
y_val = compute_y(z, r1, r2, A)
return ((y_val / C(z)) ** 1.5) * S(z) + A * np.sqrt(y_val) - np.sqrt(mu) * delta_t
# Lambert solver
def lambert(r1, r2, deltat):
# Step 1: Find the delta_theta between P1 and P2
r1_mag, r2_mag = np.linalg.norm(r1), np.linalg.norm(r2)
delta_theta = np.arccos(np.dot(r1, r2) / (r1_mag * r2_mag))
if np.cross(r1, r2)[2] < 0:
delta_theta = 2 * np.pi - delta_theta
# Step 2: Solve for the intermediate auxiliary variable z
A = np.sin(delta_theta) * np.sqrt(r1_mag * r2_mag / (1 - np.cos(delta_theta)))
def func(z_val):
return F_lambert(z_val, deltat, r1_mag, r2_mag, A)
def solver(func, z_min=-1e3, z_max=1e3, z_guess=0.01):
try:
z = brentq(lambda z: func(z), z_min, z_max)
except ValueError:
try:
z = newton(func, z_guess)
except RuntimeError:
z = least_squares(func, z_guess).x[0]
return z
z = solver(func)
y_val = compute_y(z, r1_mag, r2_mag, A)
# Step 3: Compute Lagrange coefficients f, g, gdot
f = 1 - y_val / r1_mag
g = A * np.sqrt(y_val / mu)
gdot = 1 - y_val / r2_mag
v1 = (r2 - f * r1) / g
v2 = (gdot * r2 - r1) / g
return v1, v2
v1_d, v2_a = lambert(r_earth, r_jupiter, delta_t_earth_jupiter)
v2_d, v3_a = lambert(r_jupiter, r_saturn, delta_t_jupiter_saturn)
v3_d, v4_a = lambert(r_saturn, r_uranus, delta_t_saturn_uranus)
v4_d, v5_a = lambert(r_uranus, r_neptune, delta_t_uranus_neptune)
print(v1_d, v2_a)
print(v2_d, v3_a)
print(v3_d, v4_a)
print(v4_d, v5_a)
[16.74455318 34.95333345 2.38458348] [-9.58617521 -1.22743495 -0.35734302] [-16.60866876 -10.4653879 0.87506523] [ -8.892451 -12.58488149 0.6225087 ] [ 3.38230315 -18.31306338 -0.30955557] [ 7.28171252 -14.19004446 -0.45029021] [ 13.96910355 -17.41059278 0.79967209] [ 13.94444629 -16.42954098 0.79276345]
In [22]:
print(np.linalg.norm(v1_d), np.linalg.norm(v2_a))
print(np.linalg.norm(v2_d), np.linalg.norm(v3_a))
print(np.linalg.norm(v3_d), np.linalg.norm(v4_a))
print(np.linalg.norm(v4_d), np.linalg.norm(v5_a))
38.83042387903412 9.671041604457056 19.65039341010871 15.422141356055736 18.62536146178326 15.955671725777824 22.336160600690874 21.563994831378054
4) We can now work out $\Delta V_{injection}$, the initial injection burn from $P_1$ into transfer orbit,¶
In [24]:
deltaV_injection = v1_d - v_earth
deltaV_injection = np.linalg.norm(deltaV_injection)
deltaV_injection
Out[24]:
9.777779726682427
Flyby Jupiter¶
5) Find $\Delta V_{p\,{(flyby)}}$, we must first find the hyperbolic excess velocities $\{\vec{\mathbf{v}}_{\infty_{in}},\vec{\mathbf{v}}_{\infty_{out}} \}$ at each flyby planet¶
6) Find the planet's inbound $a_{in}$ and outbound $a_{out}$ semi-major axis¶
7) Find its deflection (turning) angle $\delta$ from entering SOI to exiting SOI¶
Here, we can utilise the equation relating deflection (turning) angle $\delta$ to the outbound eccentricity $e_{out}$ of hyperbolic trajectory,¶
8) Hence, the radius of the closest approach to the planet $r_{p\,{(flyby)}}$ can be calculated¶
9) We now have all the required parameter to equate $\Delta V_{p\,{(flyby)}}$,¶
In [32]:
R_jupiter = 69911 # km
R_saturn = 58232 # km
R_uranus = 25362 # km
R_neptune = 24622 # km
R_jupiter_SOI = 48.2e6 # km
R_saturn_SOI = 54.5e6 # km
R_uranus_SOI = 51.9e6 # km
R_neptune_SOI = 86.2e6 # km
In [33]:
def find_deltaV_flyby(v_in, v_out, v_planet, mu_planet, R_planet, R_SOI):
vinf_in = v_in - v_planet
vinf_out = v_out - v_planet
a_in = -mu_planet / np.linalg.norm(vinf_in) ** 2
a_out = -mu_planet / np.linalg.norm(vinf_out) ** 2
deflection = np.arccos(
np.dot(vinf_in, vinf_out) / (np.linalg.norm(vinf_in) * np.linalg.norm(vinf_out))
)
# e_out = 1/np.sin(deflection/2)
def f(e_out, a_out, a_in, deflection):
return (a_out / a_in * (e_out - 1)) * np.sin(
deflection - np.arcsin(1 / e_out)
) - 1
def solver_e_out(a_out, a_in, deflection, e_guess=1.1, e_min=1.01, e_max=10.0):
try:
e_out = brentq(lambda e: f(e, a_out, a_in, deflection), e_min, e_max)
except ValueError:
try:
e_out = newton(lambda e: f(e, a_out, a_in, deflection), e_guess)
except RuntimeError:
e_out = least_squares(
lambda e: f(e, a_out, a_in, deflection), e_guess
).x[0]
return e_out
e_out = solver_e_out(a_out, a_in, deflection)
rp = a_out * (1 - e_out)
if rp < R_planet and rp > R_SOI:
deltaV_flyby = 1e99
else:
vp_esc = np.sqrt(2 * mu_planet / rp)
vp = vp_flyby_in = np.sqrt(np.linalg.norm(vinf_in) ** 2 + (2 * mu_planet) / rp)
if vp_flyby_in > vp_esc:
deltaV_flyby = abs(
np.sqrt(np.linalg.norm(vinf_out) ** 2 + (2 * mu_planet) / rp)
- np.sqrt(np.linalg.norm(vinf_in) ** 2 + (2 * mu_planet) / rp)
)
else:
deltaV_flyby = 1e99
return deltaV_flyby, rp, vp, vp_esc, np.rad2deg(deflection)
deltaV_jupiter, rp_jupiter, vp_jupiter, vp_esc_jupiter, deflection_jupiter = (
find_deltaV_flyby(v2_a, v2_d, v_jupiter, mu_jupiter, R_jupiter, R_jupiter_SOI)
)
display(deltaV_jupiter, rp_jupiter, vp_jupiter, vp_esc_jupiter, deflection_jupiter)
print()
deltaV_saturn, rp_saturn, vp_saturn, vp_esc_saturn, deflection_saturn = (
find_deltaV_flyby(v3_a, v3_d, v_saturn, mu_saturn, R_saturn, R_saturn_SOI)
)
display(deltaV_saturn, rp_saturn, vp_saturn, vp_esc_saturn, deflection_saturn)
print()
deltaV_uranus, rp_uranus, vp_uranus, vp_esc_uranus, deflection_uranus = (
find_deltaV_flyby(v4_a, v4_d, v_uranus, mu_uranus, R_uranus, R_uranus_SOI)
)
display(deltaV_uranus, rp_uranus, vp_uranus, vp_esc_uranus, deflection_uranus)
0.03485226340972147
2220733.0033500125
13.264161047559059
10.682600912576174
96.2960163290349
1.0288185212545642
403827.9428746656
17.457212514618856
13.706139564574247
85.73086421478894
4.086655221591039
133605.1322173773
15.891238081241008
9.313017043900006
22.164348557554597
10) Finally, to summate all the individual $\Delta V$ components of the initial injection and subsequent flybys we can equate,¶
In [35]:
deltaV_total = deltaV_injection + deltaV_jupiter + deltaV_saturn + deltaV_uranus
deltaV_total
Out[35]:
14.928105732937752
In [36]:
(
date_departure_mjd,
date_flyby_jupiter_mjd,
date_flyby_saturn_mjd,
date_flyby_uranus_mjd,
date_flyby_neptune_mjd,
)
Out[36]:
(43378.478599537164, 44063.93739583343, 44842.142326388974, 46666.74984953692, 47763.164305555634)
11) From the calculated mission $\Delta V$, we can find the fuel mass $M_f$ used.¶
In [39]:
M0_Voyager = 2066 # kg
M0_propulsion_module = 1207 # kg
Mf_propulsion_module = 1039 # kg
M0_mission_module = 825 # kg
Mf_mission_module = 100 # kg
Mf_total_NASA = Mf_propulsion_module + Mf_mission_module
In [40]:
F = 68054 # N
t = 43 # s
g0 = 9.80665 # m/s2
# The Isp of the actual Voyager I mission - propulsion module
Isp_NASA = (F * t) / (Mf_propulsion_module * g0)
Isp_NASA, "s"
Out[40]:
(287.2009612891239, 's')
Specific impulse ($I_{sp}$) is a measure of a rocket engine's efficiency, defined as the impulse (thrust per unit weight flow of propellant) delivered per unit of propellant consumed.
My Real Simulation Data¶
In [43]:
def find_Isp(Mf, M0, deltaV):
Isp = (deltaV * 1e3) / (np.log(M0 / (M0 - Mf)) * g0)
return Isp
Isp_my_sim = find_Isp(Mf_total_NASA, M0_Voyager, deltaV_total)
Isp_my_sim, "s"
Out[43]:
(1899.4417190615275, 's')
The difference in $I_{sp}$ of NASA data and my_sim data is due to the inherit assumption in my model (Method of Patched Conics and Lambert Transfers) -- also the unpowered gravitational assist assumption.
Also, due to the fact NASA used a 2 stage system with the propulsion system ejected after injection into Earth-Jupiter transfer orbit. This left a 825kg mission module with around 100kg of fuel on it.
In [45]:
Vex = Isp_my_sim * g0 / 1000 # km/s
Vex, "km/s"
Out[45]:
(18.627160134234728, 'km/s')
In [46]:
M0_Voyager = M0 = 2066 # kg
In [47]:
def compute_fuel_mass(M0, Vex, deltaV):
Mf = M0 * (1 - np.exp(-deltaV / Vex))
return Mf
fuel_mass_total = compute_fuel_mass(M0, Vex, deltaV_total)
print(
f"The mission deltaV of {deltaV_total:.4f} km/s caused {fuel_mass_total:.2f} kg of fuel used."
)
The mission deltaV of 14.9281 km/s caused 1139.00 kg of fuel used.
Function: find_deltaV_mission()¶
In [50]:
def find_deltaV_mission(
date_departure_mjd,
date_flyby_jupiter_mjd,
date_flyby_saturn_mjd,
date_flyby_uranus_mjd,
date_flyby_neptune_mjd,
):
# Step 1 - Find delta t
delta_t_earth_jupiter = find_tof(date_departure_mjd, date_flyby_jupiter_mjd)
delta_t_jupiter_saturn = find_tof(date_flyby_jupiter_mjd, date_flyby_saturn_mjd)
delta_t_saturn_uranus = find_tof(date_flyby_saturn_mjd, date_flyby_uranus_mjd)
delta_t_uranus_neptune = find_tof(date_flyby_uranus_mjd, date_flyby_neptune_mjd)
tof = (
delta_t_earth_jupiter
+ delta_t_jupiter_saturn
+ delta_t_saturn_uranus
+ delta_t_uranus_neptune
) # s
# Step 2 - Find planetary state vectors
scale = "tt"
r_earth, v_earth = get_planet_state_vector(Earth, date_departure_mjd, scale)
r_jupiter, v_jupiter = get_planet_state_vector(
Jupiter, date_flyby_jupiter_mjd, scale
)
r_saturn, v_saturn = get_planet_state_vector(Saturn, date_flyby_saturn_mjd, scale)
r_uranus, v_uranus = get_planet_state_vector(Uranus, date_flyby_uranus_mjd, scale)
r_neptune, v_neptune = get_planet_state_vector(
Neptune, date_flyby_neptune_mjd, scale
)
r_earth, v_earth = r_earth.value, v_earth.value
r_jupiter, v_jupiter = r_jupiter.value, v_jupiter.value
r_saturn, v_saturn = r_saturn.value, v_saturn.value
r_uranus, v_uranus = r_uranus.value, v_uranus.value
r_neptune, v_neptune = r_neptune.value, v_neptune.value
# Step 3 - Define trajectory fully using Lambert Transfer
v1_d, v2_a = lambert(r_earth, r_jupiter, delta_t_earth_jupiter)
v2_d, v3_a = lambert(r_jupiter, r_saturn, delta_t_jupiter_saturn)
v3_d, v4_a = lambert(r_saturn, r_uranus, delta_t_saturn_uranus)
v4_d, v5_a = lambert(r_uranus, r_neptune, delta_t_uranus_neptune)
# Step 4 - Injection deltaV
deltaV_injection = v1_d - v_earth
deltaV_injection = np.linalg.norm(deltaV_injection)
# Step 5 - Find deltaV of flyby
deltaV_jupiter, rp_jupiter, vp_jupiter, vp_esc_jupiter, deflection_jupiter = (
find_deltaV_flyby(v2_a, v2_d, v_jupiter, mu_jupiter, R_jupiter, R_jupiter_SOI)
)
deltaV_saturn, rp_saturn, vp_saturn, vp_esc_saturn, deflection_saturn = (
find_deltaV_flyby(v3_a, v3_d, v_saturn, mu_saturn, R_saturn, R_saturn_SOI)
)
deltaV_uranus, rp_uranus, vp_uranus, vp_esc_uranus, deflection_uranus = (
find_deltaV_flyby(v4_a, v4_d, v_uranus, mu_uranus, R_uranus, R_uranus_SOI)
)
# Step 6 - Find total deltaV of mission
deltaV_total = deltaV_injection + deltaV_jupiter + deltaV_saturn + deltaV_uranus
if (
delta_t_earth_jupiter <= 0
or delta_t_jupiter_saturn <= 0
or delta_t_saturn_uranus <= 0
or delta_t_uranus_neptune <= 0
):
deltaV_total = np.inf
# Step 7 - Find the fuel mass used
fuel_mass_total = compute_fuel_mass(M0, Vex, deltaV_total)
data_dict = {
"deltaV_injection (km/s)": deltaV_injection,
"deltaV_jupiter (km/s)": deltaV_jupiter,
"rp_jupiter (km)": rp_jupiter,
"vp_jupiter (km/s)": vp_jupiter,
"vp_esc_jupiter (km/s)": vp_esc_jupiter,
"deflection_jupiter (deg)": deflection_jupiter,
"deltaV_saturn (km/s)": deltaV_saturn,
"rp_saturn (km)": rp_saturn,
"vp_saturn (km/s)": vp_saturn,
"vp_esc_saturn (km/s)": vp_esc_saturn,
"deflection_saturn (deg)": deflection_saturn,
"deltaV_uranus (km/s)": deltaV_uranus,
"rp_uranus (km)": rp_uranus,
"vp_uranus (km/s)": vp_uranus,
"vp_esc_uranus (km/s)": vp_esc_uranus,
"deflection_uranus (deg)": deflection_uranus,
"tof (days)": tof / (60 * 60 * 24),
}
return deltaV_total, fuel_mass_total, data_dict
deltaV_REAL, fuel_mass_REAL, data_dict_REAL = find_deltaV_mission(
date_departure_mjd,
date_flyby_jupiter_mjd,
date_flyby_saturn_mjd,
date_flyby_uranus_mjd,
date_flyby_neptune_mjd,
)
deltaV_REAL, fuel_mass_REAL, data_dict_REAL
Out[50]:
(14.928105732937752,
1139.0,
{'deltaV_injection (km/s)': 9.777779726682427,
'deltaV_jupiter (km/s)': 0.03485226340972147,
'rp_jupiter (km)': 2220733.0033500125,
'vp_jupiter (km/s)': 13.264161047559059,
'vp_esc_jupiter (km/s)': 10.682600912576174,
'deflection_jupiter (deg)': 96.2960163290349,
'deltaV_saturn (km/s)': 1.0288185212545642,
'rp_saturn (km)': 403827.9428746656,
'vp_saturn (km/s)': 17.457212514618856,
'vp_esc_saturn (km/s)': 13.706139564574247,
'deflection_saturn (deg)': 85.73086421478894,
'deltaV_uranus (km/s)': 4.086655221591039,
'rp_uranus (km)': 133605.1322173773,
'vp_uranus (km/s)': 15.891238081241008,
'vp_esc_uranus (km/s)': 9.313017043900006,
'deflection_uranus (deg)': 22.164348557554597,
'tof (days)': 4384.68570601847})
In [51]:
dates_REAL = [
date_departure_mjd,
date_flyby_jupiter_mjd,
date_flyby_saturn_mjd,
date_flyby_uranus_mjd,
date_flyby_neptune_mjd,
]
dates_REAL
Out[51]:
[43378.478599537164, 44063.93739583343, 44842.142326388974, 46666.74984953692, 47763.164305555634]
Optimisation Algorithms¶
In [54]:
def MJD_to_TT_calander_date(T_mjd):
T_TT = Time(T_mjd, format="mjd", scale="tt") # Time object: scale='tt' format='mjd'
return T_TT.iso
MJD_to_TT_calander_date(date_departure_mjd)
Out[54]:
'1977-08-23 11:29:11.000'
In [55]:
convert_ET(MJD_to_TT_calander_date(date_departure_mjd))["UTC"]
Out[55]:
<Time object: scale='utc' format='iso' value=1977-08-23 11:28:22.816>
In [56]:
change_in_time_days = 365
Brute force Algorithm¶
1.1) Find the $\Delta t$ for each mission trajectory¶
In [59]:
print("Earth-Jupiter:", delta_t_earth_jupiter / (60 * 60 * 24), "days")
print("Jupiter-Saturn:", delta_t_jupiter_saturn / (60 * 60 * 24), "days")
print("Saturn-Uranus:", delta_t_saturn_uranus / (60 * 60 * 24), "days")
print("Uranus-Neptune:", delta_t_uranus_neptune / (60 * 60 * 24), "days")
Earth-Jupiter: 685.4587962962687 days Jupiter-Saturn: 778.2049305555411 days Saturn-Uranus: 1824.607523147948 days Uranus-Neptune: 1096.4144560187124 days
1.2) Define the range for each date allowed for algorithm¶
In [61]:
population = spacing = 15
print("No. of Iterations:", spacing**5)
No. of Iterations: 759375
In [62]:
time_range_LIST = np.linspace(-change_in_time_days, change_in_time_days, spacing)
time_range_LIST
Out[62]:
array([-365. , -312.85714286, -260.71428571, -208.57142857,
-156.42857143, -104.28571429, -52.14285714, 0. ,
52.14285714, 104.28571429, 156.42857143, 208.57142857,
260.71428571, 312.85714286, 365. ])
In [63]:
departure_time_range_LIST = np.linspace(
-change_in_time_days, change_in_time_days, spacing
)
departure_time_range_LIST
Out[63]:
array([-365. , -312.85714286, -260.71428571, -208.57142857,
-156.42857143, -104.28571429, -52.14285714, 0. ,
52.14285714, 104.28571429, 156.42857143, 208.57142857,
260.71428571, 312.85714286, 365. ])
In [64]:
flyby_jupiter_time_range_LIST = np.linspace(
-change_in_time_days, change_in_time_days, spacing
)
flyby_jupiter_time_range_LIST
Out[64]:
array([-365. , -312.85714286, -260.71428571, -208.57142857,
-156.42857143, -104.28571429, -52.14285714, 0. ,
52.14285714, 104.28571429, 156.42857143, 208.57142857,
260.71428571, 312.85714286, 365. ])
In [65]:
flyby_saturn_time_range_LIST = np.linspace(
-change_in_time_days, change_in_time_days, spacing
)
flyby_saturn_time_range_LIST
Out[65]:
array([-365. , -312.85714286, -260.71428571, -208.57142857,
-156.42857143, -104.28571429, -52.14285714, 0. ,
52.14285714, 104.28571429, 156.42857143, 208.57142857,
260.71428571, 312.85714286, 365. ])
In [66]:
flyby_uranus_time_range_LIST = np.linspace(
-change_in_time_days, change_in_time_days, spacing
)
flyby_uranus_time_range_LIST
Out[66]:
array([-365. , -312.85714286, -260.71428571, -208.57142857,
-156.42857143, -104.28571429, -52.14285714, 0. ,
52.14285714, 104.28571429, 156.42857143, 208.57142857,
260.71428571, 312.85714286, 365. ])
In [67]:
flyby_neptune_time_range_LIST = np.linspace(
-change_in_time_days, change_in_time_days, spacing
)
flyby_neptune_time_range_LIST
Out[67]:
array([-365. , -312.85714286, -260.71428571, -208.57142857,
-156.42857143, -104.28571429, -52.14285714, 0. ,
52.14285714, 104.28571429, 156.42857143, 208.57142857,
260.71428571, 312.85714286, 365. ])
In [68]:
deltaV_REAL # , fuel_mass_REAL
Out[68]:
14.928105732937752
1.3) Run Algorithm¶
In [70]:
%%time
# warnings.filterwarnings('ignore')
if global_flag:
deltaV_list_BRUTE = []
time_range_list_index = []
dates_list = []
# Start timing
start_time = time.time()
# Nested loops for all 5 mission epochs
for i in tqdm(range(0, spacing)):
date_departure_mjd_BRUTE = date_departure_mjd + time_range_LIST[i]
for j in range(0, spacing):
date_flyby_jupiter_mjd_BRUTE = date_flyby_jupiter_mjd + time_range_LIST[j]
for k in range(0, spacing):
date_flyby_saturn_mjd_BRUTE = date_flyby_saturn_mjd + time_range_LIST[k]
for l in range(0, spacing):
date_flyby_uranus_mjd_BRUTE = (
date_flyby_uranus_mjd + time_range_LIST[l]
)
for m in range(0, spacing):
date_flyby_neptune_mjd_BRUTE = (
date_flyby_neptune_mjd + time_range_LIST[m]
)
# Store index positions
time_range_list_index.append([i, j, k, l, m])
try:
# Calculate delta-V for the trajectory
deltaV_mission = find_deltaV_mission(
date_departure_mjd_BRUTE,
date_flyby_jupiter_mjd_BRUTE,
date_flyby_saturn_mjd_BRUTE,
date_flyby_uranus_mjd_BRUTE,
date_flyby_neptune_mjd_BRUTE,
)
# Append results
deltaV_list_BRUTE.append(deltaV_mission)
dates_list.append(
[
date_departure_mjd_BRUTE,
date_flyby_jupiter_mjd_BRUTE,
date_flyby_saturn_mjd_BRUTE,
date_flyby_uranus_mjd_BRUTE,
date_flyby_neptune_mjd_BRUTE,
deltaV_mission,
]
)
except Exception:
deltaV_list_BRUTE.append(
[1e99, 1e99, {"temp": 1e99}]
) # Handle failures gracefully
# End timing
end_time = time.time()
execution_time_BRUTE_FORCE = end_time - start_time
# Print execution time
print(f"\nExecution Time: {execution_time_BRUTE_FORCE:.4f} seconds")
CPU times: user 2 µs, sys: 0 ns, total: 2 µs Wall time: 3.81 µs
For spacing = For spacing = 15 --> Execution Time = 1h 15min
In [72]:
# deltaV_list_BRUTE
# Contour plot of the problem space --> look at local minima
# Holo View
# deltaV_list_BRUTE
# dates_list
In [73]:
if global_flag:
filename_brute = output_pickle.get_path("brute_results_V2.pkl")
results_brute = {
"deltaV_list_BRUTE": deltaV_list_BRUTE,
"time_range_list_index": time_range_list_index,
"dates_list": dates_list,
"execution_time_BRUTE_FORCE": execution_time_BRUTE_FORCE,
}
with open(filename_brute, "wb") as f:
pickle.dump(results_brute, f)
In [74]:
filename_brute = output_pickle.get_path("brute_results_V2.pkl")
with open(filename_brute, "rb") as f:
data = pickle.load(f)
deltaV_list_BRUTE = data["deltaV_list_BRUTE"]
time_range_list_index = data["time_range_list_index"]
dates_list = data["dates_list"]
execution_time_BRUTE_FORCE = data["execution_time_BRUTE_FORCE"]
if not global_flag:
print(f"Execution Time: {execution_time_BRUTE_FORCE:.4f} seconds")
Execution Time: 4812.2316 seconds
In [75]:
time.sleep(1) # Give the kernel a moment to catch up
In [76]:
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN = min(
deltaV_list_BRUTE, key=lambda x: x[0]
) # min(deltaV_list_BRUTE)
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN
Out[76]:
(9.965627065328311,
856.0135081245757,
{'deltaV_injection (km/s)': 9.777779726682427,
'deltaV_jupiter (km/s)': 0.03485226340972147,
'rp_jupiter (km)': 2220733.0033500125,
'vp_jupiter (km/s)': 13.264161047559059,
'vp_esc_jupiter (km/s)': 10.682600912576174,
'deflection_jupiter (deg)': 96.2960163290349,
'deltaV_saturn (km/s)': 0.0819711822100011,
'rp_saturn (km)': 385095.0255481553,
'vp_saturn (km/s)': 17.717014962912074,
'vp_esc_saturn (km/s)': 14.035548160975852,
'deflection_saturn (deg)': 85.03234833966007,
'deltaV_uranus (km/s)': 0.07102389302616174,
'rp_uranus (km)': 121875.94940717192,
'vp_uranus (km/s)': 17.652953215434174,
'vp_esc_uranus (km/s)': 9.750860929686262,
'deflection_uranus (deg)': 22.953523176610197,
'tof (days)': 4384.68570601847})
In [77]:
data_dict_BRUTE_MIN["tof (days)"]
Out[77]:
4384.68570601847
In [78]:
# Find the index of the minimum value
min_index = deltaV_list_BRUTE.index(min(deltaV_list_BRUTE, key=lambda x: x[0]))
min_index
Out[78]:
379627
In [79]:
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN = deltaV_list_BRUTE[
min_index
]
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN
Out[79]:
(9.965627065328311,
856.0135081245757,
{'deltaV_injection (km/s)': 9.777779726682427,
'deltaV_jupiter (km/s)': 0.03485226340972147,
'rp_jupiter (km)': 2220733.0033500125,
'vp_jupiter (km/s)': 13.264161047559059,
'vp_esc_jupiter (km/s)': 10.682600912576174,
'deflection_jupiter (deg)': 96.2960163290349,
'deltaV_saturn (km/s)': 0.0819711822100011,
'rp_saturn (km)': 385095.0255481553,
'vp_saturn (km/s)': 17.717014962912074,
'vp_esc_saturn (km/s)': 14.035548160975852,
'deflection_saturn (deg)': 85.03234833966007,
'deltaV_uranus (km/s)': 0.07102389302616174,
'rp_uranus (km)': 121875.94940717192,
'vp_uranus (km/s)': 17.652953215434174,
'vp_esc_uranus (km/s)': 9.750860929686262,
'deflection_uranus (deg)': 22.953523176610197,
'tof (days)': 4384.68570601847})
In [80]:
change_in_time = time_range_list_index[min_index]
change_in_time
Out[80]:
[7, 7, 7, 3, 7]
In [81]:
date_departure_mjd_BRUTE_MIN = date_departure_mjd + time_range_LIST[change_in_time[0]]
date_flyby_jupiter_mjd_BRUTE_MIN = (
date_flyby_jupiter_mjd + time_range_LIST[change_in_time[1]]
)
date_flyby_saturn_mjd_BRUTE_MIN = (
date_flyby_saturn_mjd + time_range_LIST[change_in_time[2]]
)
date_flyby_uranus_mjd_BRUTE_MIN = (
date_flyby_uranus_mjd + time_range_LIST[change_in_time[3]]
)
date_flyby_neptune_mjd_BRUTE_MIN = (
date_flyby_neptune_mjd + time_range_LIST[change_in_time[4]]
)
dates_list_BRUTE_MIN = [
date_departure_mjd_BRUTE_MIN,
date_flyby_jupiter_mjd_BRUTE_MIN,
date_flyby_saturn_mjd_BRUTE_MIN,
date_flyby_uranus_mjd_BRUTE_MIN,
date_flyby_neptune_mjd_BRUTE_MIN,
]
dates_list_BRUTE_MIN
Out[81]:
[43378.478599537164, 44063.93739583343, 44842.142326388974, 46458.178420965494, 47763.164305555634]
In [82]:
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN = find_deltaV_mission(
*dates_list_BRUTE_MIN
)
deltaV_BRUTE_MIN, fuel_mass_BRUTE_MIN, data_dict_BRUTE_MIN
Out[82]:
(9.965627065328311,
856.0135081245757,
{'deltaV_injection (km/s)': 9.777779726682427,
'deltaV_jupiter (km/s)': 0.03485226340972147,
'rp_jupiter (km)': 2220733.0033500125,
'vp_jupiter (km/s)': 13.264161047559059,
'vp_esc_jupiter (km/s)': 10.682600912576174,
'deflection_jupiter (deg)': 96.2960163290349,
'deltaV_saturn (km/s)': 0.0819711822100011,
'rp_saturn (km)': 385095.0255481553,
'vp_saturn (km/s)': 17.717014962912074,
'vp_esc_saturn (km/s)': 14.035548160975852,
'deflection_saturn (deg)': 85.03234833966007,
'deltaV_uranus (km/s)': 0.07102389302616174,
'rp_uranus (km)': 121875.94940717192,
'vp_uranus (km/s)': 17.652953215434174,
'vp_esc_uranus (km/s)': 9.750860929686262,
'deflection_uranus (deg)': 22.953523176610197,
'tof (days)': 4384.68570601847})
In [83]:
change_in_dates_list_BRUTE_MIN = [
(date_departure_mjd_BRUTE_MIN - date_departure_mjd),
(date_flyby_jupiter_mjd_BRUTE_MIN - date_flyby_jupiter_mjd),
(date_flyby_saturn_mjd_BRUTE_MIN - date_flyby_saturn_mjd),
(date_flyby_uranus_mjd_BRUTE_MIN - date_flyby_uranus_mjd),
(date_flyby_neptune_mjd_BRUTE_MIN - date_flyby_neptune_mjd),
]
change_in_dates_list_BRUTE_MIN
Out[83]:
[0.0, 0.0, 0.0, -208.57142857142753, 0.0]
In [84]:
print(
"Earth-Jupiter:",
date_flyby_jupiter_mjd_BRUTE_MIN - date_departure_mjd_BRUTE_MIN,
"days",
)
print(
"Jupiter-Saturn:",
date_flyby_saturn_mjd_BRUTE_MIN - date_flyby_jupiter_mjd_BRUTE_MIN,
"days",
)
print(
"Saturn-Uranus:",
date_flyby_uranus_mjd_BRUTE_MIN - date_flyby_saturn_mjd_BRUTE_MIN,
"days",
)
print(
"Uranus-Neptune:",
date_flyby_neptune_mjd_BRUTE_MIN - date_flyby_uranus_mjd_BRUTE_MIN,
"days",
)
Earth-Jupiter: 685.4587962962687 days Jupiter-Saturn: 778.2049305555411 days Saturn-Uranus: 1616.0360945765206 days Uranus-Neptune: 1304.98588459014 days
In [85]:
date_departure_BRUTE_MIN = MJD_to_TT_calander_date(date_departure_mjd_BRUTE_MIN)
date_flyby_jupiter_BRUTE_MIN = MJD_to_TT_calander_date(date_flyby_jupiter_mjd_BRUTE_MIN)
date_flyby_saturn_BRUTE_MIN = MJD_to_TT_calander_date(date_flyby_saturn_mjd_BRUTE_MIN)
date_flyby_uranus_BRUTE_MIN = MJD_to_TT_calander_date(
date_flyby_uranus_mjd_BRUTE_MIN
) # Fixed
date_flyby_neptune_BRUTE_MIN = MJD_to_TT_calander_date(
date_flyby_neptune_mjd_BRUTE_MIN
) # Fixed
print("Date of Departure (Min):", date_departure_BRUTE_MIN)
print("Date of Jupiter Flyby (Min):", date_flyby_jupiter_BRUTE_MIN)
print("Date of Saturn Flyby (Min):", date_flyby_saturn_BRUTE_MIN)
print("Date of Uranus Flyby (Min):", date_flyby_uranus_BRUTE_MIN) # Fixed
print("Date of Neptune Flyby (Min):", date_flyby_neptune_BRUTE_MIN) # Fixed
Date of Departure (Min): 1977-08-23 11:29:11.000 Date of Jupiter Flyby (Min): 1979-07-09 22:29:51.000 Date of Saturn Flyby (Min): 1981-08-26 03:24:57.000 Date of Uranus Flyby (Min): 1986-01-28 04:16:55.571 Date of Neptune Flyby (Min): 1989-08-25 03:56:36.000
In [86]:
# Extracting MJD values from the brute force results
departure_mjd = [dates_list[num][0] for num in range(0, len(dates_list))]
jupiter_flyby_mjd = [dates_list[num][1] for num in range(0, len(dates_list))]
saturn_flyby_mjd = [dates_list[num][2] for num in range(0, len(dates_list))]
uranus_flyby_mjd = [dates_list[num][3] for num in range(0, len(dates_list))]
neptune_flyby_mjd = [dates_list[num][4] for num in range(0, len(dates_list))]
deltaV_values = [dates_list[num][5][0] for num in range(0, len(dates_list))]
departure_mjd_change_in_time = [
(element - date_departure_mjd) for element in departure_mjd
]
jupiter_flyby_mjd_change_in_time = [
(element - date_flyby_jupiter_mjd) for element in jupiter_flyby_mjd
]
saturn_flyby_mjd_change_in_time = [
(element - date_flyby_saturn_mjd) for element in saturn_flyby_mjd
]
uranus_flyby_mjd_change_in_time = [
(element - date_flyby_uranus_mjd) for element in uranus_flyby_mjd
]
neptune_flyby_mjd_change_in_time = [
(element - date_flyby_neptune_mjd) for element in neptune_flyby_mjd
]
1.4) Contour Graph of Problem Space¶
In [88]:
# Create a grid for the contour plot
max_deltav_threshold = 15 # km/s (adjust as needed)
valid_indices = [
i for i in range(len(deltaV_values)) if deltaV_values[i] < max_deltav_threshold
]
# Extract valid data
valid_jupiter_mjd = [jupiter_flyby_mjd_change_in_time[i] for i in valid_indices]
valid_saturn_mjd = [saturn_flyby_mjd_change_in_time[i] for i in valid_indices]
valid_deltaV = [deltaV_values[i] for i in valid_indices]
# Grid sizes
jupiter_range = np.linspace(min(valid_jupiter_mjd), max(valid_jupiter_mjd), 1000)
saturn_range = np.linspace(min(valid_saturn_mjd), max(valid_saturn_mjd), 1000)
# Create meshgrid
X, Y = np.meshgrid(jupiter_range, saturn_range)
# Interpolate deltaV values to the grid using linear interpolation
Z = griddata(
(valid_jupiter_mjd, valid_saturn_mjd), valid_deltaV, (X, Y), method="linear"
)
# Create the contour plot
fig, ax = plt.subplots(figsize=(10, 8))
# LaTeX-style text and font for better readability
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
# Set up a colormap that emphasises the minimum Delta-V
cmap = plt.cm.viridis_r
norm = colors.Normalize(vmin=np.nanmin(Z), vmax=np.nanmax(Z))
# Plot the contour
# Define contour levels within the specified range
levels = np.linspace(9, max_deltav_threshold, 20) # 20 levels from 9 to max threshold
# Filled contour plot
contour = ax.contourf(
X,
Y,
Z,
levels=levels,
cmap=cmap,
norm=colors.Normalize(vmin=9, vmax=max_deltav_threshold),
)
cbar = plt.colorbar(contour, ax=ax, label=r"$\Delta V$ (km/s)")
cbar.ax.tick_params(labelsize=12)
contour_lines = ax.contour(
X, Y, Z, levels=10, colors="white", linewidths=0.5, alpha=0.7
)
ax.clabel(contour_lines, inline=True, fontsize=8, fmt=r"%.2f")
# Plot key mission points
ax.scatter(
0,
0,
color="red",
s=100,
marker="o",
edgecolors="black",
linewidth=1.2,
label=r"Real $\Delta V$: {:.2f} km/s".format(deltaV_REAL),
)
ax.scatter(
date_flyby_jupiter_mjd_BRUTE_MIN - date_flyby_jupiter_mjd,
date_flyby_saturn_mjd_BRUTE_MIN - date_flyby_saturn_mjd,
color="fuchsia",
s=200,
marker="*",
edgecolors="black",
linewidth=1.2,
label=r"Minimum $\Delta V$: {:.2f} km/s".format(deltaV_BRUTE_MIN),
)
# Set axis labels
ax.set_xlabel(r"Change in Jupiter Flyby Epoch (days)", labelpad=10)
ax.set_ylabel(r"Change in Saturn Flyby Epoch (days)", labelpad=10)
# ax.set_title(r"$\Delta V$ Contour for Jupiter \& Saturn Flybys ($\pm$365 days)", pad=15)
# Improve grid and ticks
# ax.grid(which='major', color='#DDDDDD', linewidth=0.8)
# ax.grid(which='minor', color='#EEEEEE', linestyle=':', linewidth=0.5)
ax.minorticks_on()
# Customise tick parameters
ax.tick_params(axis="both", which="major", labelsize=12)
ax.tick_params(axis="both", which="minor", labelsize=10)
plt.legend(loc="upper left")
# Adjust layout for better spacing
plt.tight_layout()
output_figures.save_file(
"voyager_II_brute_force_contour_2D.png", dpi=500, bbox_inches="tight"
)
plt.show()
plt.close()
Figure saved to: Voyager 2 - Figures/voyager_II_brute_force_contour_2D.png
In [89]:
# Filter out any extreme values
max_deltav_threshold = 20 # km/s (adjust as needed)
valid_indices = [
i for i in range(len(deltaV_values)) if deltaV_values[i] < max_deltav_threshold
]
# Extract valid data
valid_departure_mjd = [departure_mjd_change_in_time[i] for i in valid_indices]
valid_jupiter_mjd = [jupiter_flyby_mjd_change_in_time[i] for i in valid_indices]
valid_saturn_mjd = [saturn_flyby_mjd_change_in_time[i] for i in valid_indices]
valid_deltaV = [deltaV_values[i] for i in valid_indices]
# Create a 3D plot
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection="3d")
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
# Create scatter plot with color based on delta-V 'viridis_r' 'plasma_r'
scatter = ax.scatter(
valid_departure_mjd,
valid_jupiter_mjd,
valid_saturn_mjd,
c=valid_deltaV,
cmap="viridis_r",
s=5,
alpha=0.25,
norm=colors.Normalize(vmin=min(valid_deltaV), vmax=min(valid_deltaV) * 1.5),
)
# Highlight the delta-V points
ax.scatter(
[0],
[0],
[0],
color="red",
s=100,
marker="o",
edgecolors="black",
linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f"Real Delta-V: {deltaV_REAL:.2f} km/s\nReal Fuel Mass: {fuel_mass_REAL:.2f} kg",
)
ax.scatter(
[date_departure_mjd_BRUTE_MIN - date_departure_mjd],
[date_flyby_jupiter_mjd_BRUTE_MIN - date_flyby_jupiter_mjd],
[date_flyby_saturn_mjd_BRUTE_MIN - date_flyby_saturn_mjd],
color="fuchsia",
s=200,
marker="*",
edgecolors="black",
linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f"Minimum Delta-V: {deltaV_BRUTE_MIN:.2f} km/s\nMinimum Fuel Mass: {fuel_mass_BRUTE_MIN:.2f} kg",
)
date_flyby_saturn_mjd_BRUTE_MIN
# Add colorbar
cbar = plt.colorbar(scatter, ax=ax, label="Delta-V (km/s)")
# Set labels and title
ax.set_xlabel("Change in Earth Departure Epoch (days)")
ax.set_ylabel("Change in Jupiter Flyby Epoch (days)")
ax.set_zlabel("Change in Saturn Flyby Epoch (days)")
# ax.set_title("3D Visualisation of Brute Force Problem Space ($\pm$ 365 days bounds)")
# Set axis limits
limit = 375
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_zlim(-limit, limit)
# Add a legend
ax.legend(loc="best")
# Show both plots
plt.tight_layout()
output_figures.save_file(
"voyager_II_brute_force_contour_3D.png", dpi=500, bbox_inches="tight"
)
# ax.view_init(elev=30, azim=45) # Adjust elevation & azimuth for a better view
plt.show()
plt.close()
Figure saved to: Voyager 2 - Figures/voyager_II_brute_force_contour_3D.png
In [90]:
gc.collect()
Out[90]:
255
In [92]:
dates_REAL
Out[92]:
[43378.478599537164, 44063.93739583343, 44842.142326388974, 46666.74984953692, 47763.164305555634]
In [93]:
%%time
if global_flag:
iterations = []
iteration_count = []
# Initial guess for mission dates
initial_guess = np.array([dates_REAL])
def objective_function(dates):
# Unpack date variables
(
date_departure_mjd,
date_flyby_jupiter_mjd,
date_flyby_saturn_mjd,
date_flyby_uranus_mjd,
date_flyby_neptune_mjd,
) = dates
iteration_count.append(dates.copy())
# Define search bounds
lower_bound = initial_guess - 365
upper_bound = initial_guess + 365
# Check if any date is out of bounds
if np.any(dates < lower_bound) or np.any(dates > upper_bound):
return np.inf # Assign a high penalty for out-of-bounds solutions
try:
# Compute Delta-V for the given dates
deltaV, fuel_mass, data_dict = find_deltaV_mission(
date_departure_mjd,
date_flyby_jupiter_mjd,
date_flyby_saturn_mjd,
date_flyby_uranus_mjd,
date_flyby_neptune_mjd,
)
# Store the iteration history
iterations.append(dates.copy())
except Exception as e:
print(f"Error in simulation: {e}")
return np.inf # Assign a high penalty for failed evaluations
return deltaV # Minimize Delta-V
# Start timing the optimisation
start_time = time.time()
# Optimise using `fmin` (Simplex method)
optimal_dates = fmin(objective_function, initial_guess, disp=True, maxiter=2000)
iteration_count = len(iteration_count)
# Extract the optimised departure and flyby dates
(
date_departure_mjd_FMIN,
date_flyby_jupiter_mjd_FMIN,
date_flyby_saturn_mjd_FMIN,
date_flyby_uranus_mjd_FMIN,
date_flyby_neptune_mjd_FMIN,
) = optimal_dates
# Compute Delta-V at optimissed dates
deltaV_FMIN, fuel_mass_FMIN, data_dict_FMIN = find_deltaV_mission(*optimal_dates)
# Display results
print(f"\nMinimum Delta-V: {deltaV_FMIN:.4f} km/s")
print(f"Optimal Departure Date (MJD): {date_departure_mjd_FMIN:.2f}")
print(f"Optimal Jupiter Flyby Date (MJD): {date_flyby_jupiter_mjd_FMIN:.2f}")
print(f"Optimal Saturn Flyby Date (MJD): {date_flyby_saturn_mjd_FMIN:.2f}")
print(f"Optimal Uranus Flyby Date (MJD): {date_flyby_uranus_mjd_FMIN:.2f}")
print(f"Optimal Neptune Flyby Date (MJD): {date_flyby_neptune_mjd_FMIN:.2f}")
# Convert MJD dates to calendar dates
print("\nIn calendar dates:")
print(f"Departure: {MJD_to_TT_calander_date(date_departure_mjd_FMIN)}")
print(f"Jupiter flyby: {MJD_to_TT_calander_date(date_flyby_jupiter_mjd_FMIN)}")
print(f"Saturn flyby: {MJD_to_TT_calander_date(date_flyby_saturn_mjd_FMIN)}")
print(f"Uranus flyby: {MJD_to_TT_calander_date(date_flyby_uranus_mjd_FMIN)}")
print(f"Neptune flyby: {MJD_to_TT_calander_date(date_flyby_neptune_mjd_FMIN)}")
# End timing
execution_time_FMIN = time.time() - start_time
print(f"\nExecution Time: {execution_time_FMIN:.4f} seconds")
print(f"No. of Iterations: {iteration_count}")
CPU times: user 2 µs, sys: 0 ns, total: 2 µs Wall time: 4.05 µs
In [94]:
if global_flag:
filename_brute = output_pickle.get_path("fmin_results_V2.pkl")
results_brute = {
"optimal_dates": optimal_dates,
"execution_time_FMIN": execution_time_FMIN,
"iterations": iterations,
"iteration_count": iteration_count,
}
with open(filename_brute, "wb") as f:
pickle.dump(results_brute, f)
In [95]:
filename_brute = output_pickle.get_path("fmin_results_V2.pkl")
with open(filename_brute, "rb") as f:
data = pickle.load(f)
optimal_dates = data["optimal_dates"]
execution_time_FMIN = data["execution_time_FMIN"]
iterations = data["iterations"]
iteration_count = data["iteration_count"]
if not global_flag:
# Extract the optimised departure and flyby dates
(
date_departure_mjd_FMIN,
date_flyby_jupiter_mjd_FMIN,
date_flyby_saturn_mjd_FMIN,
date_flyby_uranus_mjd_FMIN,
date_flyby_neptune_mjd_FMIN,
) = optimal_dates
# Compute Delta-V at optimissed dates
deltaV_FMIN, fuel_mass_FMIN, data_dict_FMIN = find_deltaV_mission(*optimal_dates)
# Display results
print(f"\nMinimum Delta-V: {deltaV_FMIN:.4f} km/s")
print(f"Optimal Departure Date (MJD): {date_departure_mjd_FMIN:.2f}")
print(f"Optimal Jupiter Flyby Date (MJD): {date_flyby_jupiter_mjd_FMIN:.2f}")
print(f"Optimal Saturn Flyby Date (MJD): {date_flyby_saturn_mjd_FMIN:.2f}")
print(f"Optimal Uranus Flyby Date (MJD): {date_flyby_uranus_mjd_FMIN:.2f}")
print(f"Optimal Neptune Flyby Date (MJD): {date_flyby_neptune_mjd_FMIN:.2f}")
# Convert MJD dates to calendar dates
print("\nIn calendar dates:")
print(f"Departure: {MJD_to_TT_calander_date(date_departure_mjd_FMIN)}")
print(f"Jupiter flyby: {MJD_to_TT_calander_date(date_flyby_jupiter_mjd_FMIN)}")
print(f"Saturn flyby: {MJD_to_TT_calander_date(date_flyby_saturn_mjd_FMIN)}")
print(f"Uranus flyby: {MJD_to_TT_calander_date(date_flyby_uranus_mjd_FMIN)}")
print(f"Neptune flyby: {MJD_to_TT_calander_date(date_flyby_neptune_mjd_FMIN)}")
print(f"\nExecution Time: {execution_time_FMIN:.4f} seconds")
print(f"No. of Iterations: {iteration_count}")
Minimum Delta-V: 9.3713 km/s Optimal Departure Date (MJD): 43388.97 Optimal Jupiter Flyby Date (MJD): 44119.21 Optimal Saturn Flyby Date (MJD): 44971.48 Optimal Uranus Flyby Date (MJD): 46706.76 Optimal Neptune Flyby Date (MJD): 48128.16 In calendar dates: Departure: 1977-09-02 23:16:50.675 Jupiter flyby: 1979-09-03 05:06:38.181 Saturn flyby: 1982-01-02 11:28:33.046 Uranus flyby: 1986-10-03 18:21:11.480 Neptune flyby: 1990-08-25 03:56:35.928 Execution Time: 5.9432 seconds No. of Iterations: 1443
In [96]:
# iterations
In [97]:
deltaV_FMIN, fuel_mass_FMIN, data_dict_FMIN = find_deltaV_mission(
date_departure_mjd_FMIN,
date_flyby_jupiter_mjd_FMIN,
date_flyby_saturn_mjd_FMIN,
date_flyby_uranus_mjd_FMIN,
date_flyby_neptune_mjd_FMIN,
)
deltaV_FMIN, fuel_mass_FMIN, data_dict_FMIN
Out[97]:
(9.37125750134365,
816.7817510508315,
{'deltaV_injection (km/s)': 9.35789550243903,
'deltaV_jupiter (km/s)': 0.0009168183205279234,
'rp_jupiter (km)': 2802517.371817503,
'vp_jupiter (km/s)': 11.834271643601488,
'vp_esc_jupiter (km/s)': 9.50935412829766,
'deflection_jupiter (deg)': 94.14765367850407,
'deltaV_saturn (km/s)': 2.6201263381153694e-11,
'rp_saturn (km)': 516702.31495290494,
'vp_saturn (km/s)': 15.304134784770053,
'vp_esc_saturn (km/s)': 12.11694480065466,
'deflection_saturn (deg)': 84.29926641857209,
'deltaV_uranus (km/s)': 0.012445180557890012,
'rp_uranus (km)': 171918.01181510068,
'vp_uranus (km/s)': 15.542222697049406,
'vp_esc_uranus (km/s)': 8.209964483111364,
'deflection_uranus (deg)': 20.471397443524044,
'tof (days)': 4739.194273761466})
In [98]:
print("Earth-Jupiter:", date_flyby_jupiter_mjd_FMIN - date_departure_mjd_FMIN, "days")
print(
"Jupiter-Saturn:", date_flyby_saturn_mjd_FMIN - date_flyby_jupiter_mjd_FMIN, "days"
)
print("Saturn-Uranus:", date_flyby_uranus_mjd_FMIN - date_flyby_saturn_mjd_FMIN, "days")
print(
"Uranus-Neptune:", date_flyby_neptune_mjd_FMIN - date_flyby_uranus_mjd_FMIN, "days"
)
Earth-Jupiter: 730.2429109525983 days Jupiter-Saturn: 852.2652183495811 days Saturn-Uranus: 1735.2865559484926 days Uranus-Neptune: 1421.3995885107943 days
In [99]:
dates_list_FMIN = [
date_departure_mjd_FMIN,
date_flyby_jupiter_mjd_FMIN,
date_flyby_saturn_mjd_FMIN,
date_flyby_uranus_mjd_FMIN,
date_flyby_neptune_mjd_FMIN,
]
dates_list_FMIN
Out[99]:
[43388.97003095798, 44119.21294191058, 44971.47816026016, 46706.76471620865, 48128.164304719445]
In [100]:
change_in_dates_list_FMIN = [
(date_departure_mjd_FMIN - date_departure_mjd),
(date_flyby_jupiter_mjd_FMIN - date_flyby_jupiter_mjd),
(date_flyby_saturn_mjd_FMIN - date_flyby_saturn_mjd),
(date_flyby_uranus_mjd_FMIN - date_flyby_uranus_mjd),
(date_flyby_neptune_mjd_FMIN - date_flyby_neptune_mjd),
]
change_in_dates_list_FMIN
Out[100]:
[10.491431420814479, 55.2755460771441, 129.33583387118415, 40.01486667172867, 364.99999916381057]
In [101]:
date_departure_FMIN = MJD_to_TT_calander_date(date_departure_mjd_FMIN)
date_flyby_jupiter_FMIN = MJD_to_TT_calander_date(date_flyby_jupiter_mjd_FMIN)
date_flyby_saturn_FMIN = MJD_to_TT_calander_date(date_flyby_saturn_mjd_FMIN)
date_flyby_uranus_FMIN = MJD_to_TT_calander_date(date_flyby_uranus_mjd_FMIN)
date_flyby_neptune_FMIN = MJD_to_TT_calander_date(date_flyby_neptune_mjd_FMIN)
print("Date of Departure (Min):", date_departure_FMIN)
print("Date of Jupiter Flyby (Min):", date_flyby_jupiter_FMIN)
print("Date of Saturn Flyby (Min):", date_flyby_saturn_FMIN)
print("Date of Uranus Flyby (Min):", date_flyby_uranus_FMIN) # Fixed
print("Date of Neptune Flyby (Min):", date_flyby_neptune_FMIN) # Fixed
Date of Departure (Min): 1977-09-02 23:16:50.675 Date of Jupiter Flyby (Min): 1979-09-03 05:06:38.181 Date of Saturn Flyby (Min): 1982-01-02 11:28:33.046 Date of Uranus Flyby (Min): 1986-10-03 18:21:11.480 Date of Neptune Flyby (Min): 1990-08-25 03:56:35.928
Contour Graph of Problem Space¶
In [103]:
# Extract the dates from the iterations list
departure_dates_FMIN = np.array([iter[0] for iter in iterations])
jupiter_flyby_dates_FMIN = np.array([iter[1] for iter in iterations])
saturn_flyby_dates_FMIN = np.array([iter[2] for iter in iterations])
# Calculate the corresponding deltaV values for the iterations
deltaV_values_FMIN = []
for dates in iterations:
try:
deltaV, _, _ = find_deltaV_mission(*dates)
except Exception:
deltaV_values_FMIN.append(np.inf)
continue
deltaV_values_FMIN.append(deltaV)
# Create a 3D plot
fig = plt.figure(figsize=(10, 8))
ax = fig.add_subplot(111, projection="3d")
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
# Scatter plot with color based on deltaV
scatter = ax.scatter(
departure_dates_FMIN - date_departure_mjd,
jupiter_flyby_dates_FMIN - date_flyby_jupiter_mjd,
saturn_flyby_dates_FMIN - date_flyby_saturn_mjd,
c=deltaV_values_FMIN,
cmap="viridis_r",
s=10,
alpha=0.7,
norm=plt.Normalize(
vmin=min(deltaV_values_FMIN), vmax=min(deltaV_values_FMIN) * 1.5
),
)
# ax.plot(departure_dates_FMIN - date_departure_mjd, jupiter_flyby_dates_FMIN - date_flyby_jupiter_mjd, saturn_flyby_dates_FMIN - date_flyby_saturn_mjd,
# linestyle='-', color='g', markersize=6, linewidth=0.5, alpha=0.5)
ax.scatter(
[0],
[0],
[0],
color="red",
s=100,
marker="o",
edgecolors="black",
linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f"Real Delta-V: {deltaV_REAL:.2f} km/s\nReal Fuel Mass: {fuel_mass_REAL:.2f} kg",
)
ax.scatter(
date_departure_mjd_FMIN - date_departure_mjd,
date_flyby_jupiter_mjd_FMIN - date_flyby_jupiter_mjd,
date_flyby_saturn_mjd_FMIN - date_flyby_saturn_mjd,
color="fuchsia",
s=200,
marker="*",
edgecolors="black",
linewidth=1.2,
zorder=3, # Higher zorder brings it to the front
label=f"Minimum Delta-V: {deltaV_FMIN:.4f} km/s\nMinimum Fuel Mass: {fuel_mass_FMIN:.2f} kg",
)
# Add colorbar
cbar = plt.colorbar(scatter, ax=ax, label="Delta-V (km/s)")
# Set labels and title
ax.set_xlabel("Change in Earth Departure Epoch (days)")
ax.set_ylabel("Change in Jupiter Flyby Epoch (days)")
ax.set_zlabel("Change in Saturn Flyby Epoch (days)")
# ax.set_title("3D Visualisation of FMIN Problem Space ($\pm$ 365 days bounds)")
ax.legend(loc="best")
# Set axis limits
# limit = 375
# ax.set_xlim(-limit, limit)
# ax.set_ylim(-limit, limit)
# ax.set_zlim(-limit, limit)
# Show the plot
plt.tight_layout()
output_figures.save_file("voyager_II_fmin_contour_3D.png", dpi=500, bbox_inches="tight")
plt.show()
plt.close()
Figure saved to: Voyager 2 - Figures/voyager_II_fmin_contour_3D.png
In [104]:
df = pd.DataFrame(
iterations,
columns=[
"Departure Date",
"Jupiter Flyby",
"Saturn Flyby",
"Uranus Flyby",
"Neptune Flyby",
],
)
# Calculate the corresponding deltaV values for the iterations
deltaV_values_FMIN = []
for dates in iterations:
try:
deltaV, _, _ = find_deltaV_mission(
*dates
) # Assuming find_deltaV_mission function is defined
except Exception:
deltaV_values_FMIN.append(np.inf)
continue
deltaV_values_FMIN.append(deltaV)
# Create an index for the iterations and add deltaV values to the DataFrame
df["deltaV"] = deltaV_values_FMIN
df["Iteration"] = df.index
# Set up the colour map for the lines based on deltaV values
norm = Normalize(vmin=np.min(deltaV_values_FMIN), vmax=np.max(deltaV_values_FMIN))
colors = plt.cm.viridis(norm(deltaV_values_FMIN))
# Create the plot
fig, ax = plt.subplots(figsize=(10, 6))
# Plotting the parallel coordinates plot with dashed lines
for i, row in df.iterrows():
ax.plot(
df.columns[:-2], row[:-2], linestyle="--", linewidth=0.5, color=colors[i]
) # Use color based on deltaV
ax.plot(
df.columns[:-2],
dates_REAL,
linestyle="--",
color="red",
linewidth=2,
label="Real Mission",
) # Highlighted line
# Highlight the last iteration
last_iteration = df.iloc[-1]
ax.plot(
df.columns[:-2],
last_iteration[:-2],
linestyle="--",
color="orange",
linewidth=2,
label="Optimised Solution (Last Iteration)",
) # Highlighted line
# Add colour bar for deltaV
sm = plt.cm.ScalarMappable(cmap="viridis_r", norm=norm)
sm.set_array([])
fig.colorbar(sm, ax=ax, label="deltaV (m/s)")
# Set title and labels
# ax.set_title("Parallel Coordinates Plot of Optimisation Over the 5 Epochs (FMIN)")
ax.set_xlabel("Epochs")
ax.set_ylabel("Dates in MJD (days)")
# Improve grid and ticks
ax.grid(which="major", color="#DDDDDD", linewidth=0.8)
ax.grid(which="minor", color="#EEEEEE", linestyle=":", linewidth=0.5)
ax.minorticks_on()
# Customise tick parameters
ax.tick_params(axis="both", which="major", labelsize=12)
ax.tick_params(axis="both", which="minor", labelsize=10)
ax.tick_params(axis="x", which="minor", length=0)
# Add legend for last iteration
ax.legend()
# Adjust layout for tightness
plt.tight_layout()
# Show the plot
plt.show()
plt.close()
In [105]:
df_chosen = df.head(50)
# Select relevant columns
df_pairwise = df_chosen[
["Departure Date", "Jupiter Flyby", "Saturn Flyby", "Uranus Flyby", "Neptune Flyby"]
]
# Compute pairwise differences
pairwise_changes = df_pairwise.diff(axis=1)
# Rename columns to match desired format
pairwise_changes.columns = [
"Departure",
"Jupiter Flyby - Departure",
"Saturn Flyby - Jupiter Flyby",
"Uranus Flyby - Saturn Flyby",
"Neptune Flyby - Uranus Flyby",
]
pairwise_changes = pairwise_changes.iloc[:, 1:]
# Compute the absolute values of the differences
# pairwise_changes = pairwise_changes.abs()
def format_annotation(val):
if val < 1e-3:
return f"{val:.2e}"
else:
return f"{int(val)}"
# Create a heatmap of pairwise changes for the last 50 iterations
# Plot the data
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
fig, ax = plt.subplots(figsize=(7, 10))
sns.heatmap(
pairwise_changes,
annot=True,
fmt=".0f",
cmap="viridis",
cbar_kws={"label": "Pairwise Change"},
annot_kws={"size": 10},
)
# Update annotations with custom format
for text in plt.gca().texts:
text.set_text(format_annotation(float(text.get_text())))
plt.xticks(rotation=20)
plt.yticks(rotation=0)
# plt.title("Heatmap of Pairwise Changes Between Flyby Epochs")
plt.xlabel("Flyby Epochs")
plt.ylabel("Iterations")
# ax.minorticks_on()
# ax.tick_params(axis='both', which='major', labelsize=12)
# ax.tick_params(axis='both', which='minor', labelsize=10)
ax.tick_params(axis="x", which="minor", length=0)
# Show the plot
plt.tight_layout()
plt.show()
plt.close()
In [106]:
df.columns[:]
Out[106]:
Index(['Departure Date', 'Jupiter Flyby', 'Saturn Flyby', 'Uranus Flyby',
'Neptune Flyby', 'deltaV', 'Iteration'],
dtype='object')
In [107]:
[
(date_departure_mjd - date_departure_mjd_FMIN),
(date_flyby_jupiter_mjd - date_flyby_jupiter_mjd_FMIN),
(date_flyby_saturn_mjd - date_flyby_saturn_mjd_FMIN),
(date_flyby_uranus_mjd - date_flyby_uranus_mjd_FMIN),
(date_flyby_neptune_mjd - date_flyby_neptune_mjd_FMIN),
]
Out[107]:
[-10.491431420814479, -55.2755460771441, -129.33583387118415, -40.01486667172867, -364.99999916381057]
Genetic Algorithm (GA)¶
In [109]:
change_in_time_days = 365
min_departure_date = date_departure_mjd - change_in_time_days
max_departure_date = date_departure_mjd + change_in_time_days
min_flyby_jupiter = date_flyby_jupiter_mjd - change_in_time_days
max_flyby_jupiter = date_flyby_jupiter_mjd + change_in_time_days
min_flyby_saturn = date_flyby_saturn_mjd - change_in_time_days
max_flyby_saturn = date_flyby_saturn_mjd + change_in_time_days
min_flyby_uranus = date_flyby_uranus_mjd - change_in_time_days
max_flyby_uranus = date_flyby_uranus_mjd + change_in_time_days
min_flyby_neptune = date_flyby_neptune_mjd - change_in_time_days
max_flyby_neptune = date_flyby_saturn_mjd + change_in_time_days
Step 1: Generate Initial Population¶
In [111]:
# Step 1: Generate Initial Population
def initial_population(n_population):
# Generates initial population within mission constraints.
population = []
population.append(dates_REAL)
# Add interval to generate population
for num in range(n_population - 1):
date_departure_pop = random.uniform(min_departure_date, max_departure_date)
date_flyby_jupiter_pop = random.uniform(min_flyby_jupiter, max_flyby_jupiter)
date_flyby_saturn_pop = random.uniform(min_flyby_saturn, max_flyby_saturn)
date_flyby_uranus_pop = random.uniform(min_flyby_uranus, max_flyby_uranus)
date_flyby_neptune_pop = random.uniform(min_flyby_neptune, max_flyby_neptune)
# Ensure chronological order of mission events
dates_ordered = np.sort(
[
date_departure_pop,
date_flyby_jupiter_pop,
date_flyby_saturn_pop,
date_flyby_uranus_pop,
date_flyby_neptune_pop,
]
)
population.append(dates_ordered)
return population
initial_population(5)
Out[111]:
[[43378.478599537164,
44063.93739583343,
44842.142326388974,
46666.74984953692,
47763.164305555634],
array([43021.39085934, 43933.77025936, 45147.45580997, 46900.77361275,
47230.68614329]),
array([43405.41902688, 44206.29392034, 45048.1038886 , 46515.75210546,
47390.21567878]),
array([43650.88322171, 43816.8355646 , 44732.24116711, 45211.63871084,
46709.07066809]),
array([43624.77203853, 43765.12390283, 44644.31815834, 46647.64340839,
46921.33900256])]
Step 2: Evaluation [Fitness Function (Objective Function)]¶
In [113]:
# Step 2: Evaluation - Fitness Function (Objective Function)
def fitness(individual):
# Evaluates fitness based on deltaV_total minimisation
try:
(
date_departure_mjd,
date_flyby_jupiter_mjd,
date_flyby_saturn_mjd,
date_flyby_uranus_mjd,
date_flyby_neptune_mjd,
) = individual
deltaV, fuel_mass, data_dict = find_deltaV_mission(*individual)
except Exception:
deltaV = 1e6 # Assign a large penalty for infeasible solutions
return 1 / (deltaV + 1e-6) # Avoid division by zero
Step 3: Sort Population by Fitness¶
In [115]:
def sort_population_by_fitness(population, fitness_scores, n_population):
sorted_indices = np.argsort(fitness_scores)[::-1] # Sort by highest fitness
population = [
population[i] for i in sorted_indices[:n_population]
] # Ensures length population is 'n_population'
return population
Step 4: Elitism (Optional)¶
In [117]:
def elitism(population, n_elite):
elite_individuals = population[:n_elite]
return elite_individuals
Step 5: Selection¶
Tournament Selection¶
In [120]:
def tournament_selection(population, fitness_scores, tournament_size=3):
# Select random individuals for the tournament
selected_indices = np.random.choice(len(population), tournament_size, replace=False)
# Find the best individual in the tournament
tournament_fitness = [fitness_scores[i] for i in selected_indices]
winner_idx = selected_indices[np.argmax(tournament_fitness)]
return population[winner_idx]
Roulette Wheel Selection¶
In [122]:
def roulette_wheel_selection(population, fitness_values):
# Calculate the normalised selection probability for each individual
selection_prob = fitness_values / sum(fitness_values)
# Calculate cumulative sum of fitness probabilities
cum_sum = selection_prob.cumsum()
# Generate a random number between 0 and 1
r = random.random()
# Find the first index where cumulative sum exceeds random number
selected_index = np.where(cum_sum >= r)[0][0]
return population[selected_index]
pop = initial_population(20)
# Evaluate fitness function
fitness_scores = np.array([fitness(ind) for ind in pop])
# Sorts population by fitness value
pop = sort_population_by_fitness(pop, fitness_scores, 20)
roulette_wheel_selection(pop, fitness_scores)
Out[122]:
array([43459.14650599, 43928.51892308, 44500.12200123, 45563.95844498,
46990.2409105 ])
Step 6: Crossover¶
Simulated Binary Crossover (SBX)¶
In [125]:
def sbx_crossover(parent_1, parent_2, eta):
offspring_1 = parent_1[:]
offspring_2 = parent_2[:]
for i in range(len(parent_1)):
u = random.random()
if u <= 0.5:
beta = (2 * u) ** (1 / (eta + 1))
else:
beta = (1 / (2 * (1 - u))) ** (1 / (eta + 1))
# Generate offspring
offspring_1[i] = 0.5 * ((1 + beta) * parent_1[i] + (1 - beta) * parent_2[i])
offspring_2[i] = 0.5 * ((1 - beta) * parent_1[i] + (1 + beta) * parent_2[i])
# Ensure chronological order (if values represent time-based sequences)
# offspring_1 = np.sort(offspring_1)
# offspring_2 = np.sort(offspring_2)
return offspring_1, offspring_2
Step 7: Mutation (Random Perturbations)¶
In [127]:
def gaussian_mutation(individual):
# Create a copy to avoid modifying the original
mutated = individual.copy()
temp = mutated
# Define bounds for each parameter
bounds = [
[min_departure_date, max_departure_date],
[min_flyby_jupiter, max_flyby_jupiter],
[min_flyby_saturn, max_flyby_saturn],
[min_flyby_uranus, max_flyby_uranus],
[min_flyby_neptune, max_flyby_neptune],
]
# Select random index to mutate
index = random.randint(0, 4)
# Calculate parameter range
param_range = bounds[index][1] - bounds[index][0]
# Apply Gaussian mutation with standard deviation proportional to parameter range
sigma = 0.05 * param_range # 5% of the parameter range
delta = random.gauss(0, sigma)
mutated[index] += delta
# Ensure value stays within bounds
if mutated[index] < bounds[index][0] or mutated[index] > bounds[index][1]:
mutated[index] = temp[index]
# Ensure sequential order
# mutated = np.sort(mutated)
return mutated
Step 8: Running GA¶
- Initialise population
- Sort population
- Start Loop
- Evaluate fitness
- Select parents (roulette wheel)
- Crossover (SBX or copy)
- Mutation
- Combine Offspring & Elite Individuals
- Filter for constraints
- Sort & truncate population
- Elitism (store best individuals)
- Does it satisfy termination condition
- If NO, then repeat loop. If YES, end loop
- Store best (optimal) individual solution
In [129]:
def run_genetic_algorithm(
n_population,
n_elite,
crossover_per,
mutation_per,
n_generations,
population_per_generation,
):
population = initial_population(n_population) # Generate initial population
# Evaluate fitness function
fitness_scores_temp = np.array([fitness(ind) for ind in population])
# Sorts population by fitness value
population_temp = sort_population_by_fitness(
population, fitness_scores_temp, n_population
)
population_per_generation.append(population_temp)
elite_individuals = []
for generation in range(0, n_generations):
start_time = time.time()
# Evaluate fitness function
fitness_scores = np.array([fitness(ind) for ind in population])
# Sorts population by fitness value
population = sort_population_by_fitness(
population, fitness_scores, n_population
)
# Print progress
best_individual = population[0]
best_deltaV, best_fuel_mass, _ = find_deltaV_mission(*best_individual)
print(
f"Generation {generation}: Best Delta-V = {best_deltaV:.6f} km/s, Fitness = {fitness(best_individual):.6f}"
)
# Apply elitism - preserve the best individuals
if n_elite != 0:
elite_individuals = elitism(population, n_elite)
# Select parents for potential crossover - select enough for the whole population
parents_list = []
for i in range(n_population):
parents_list.append(tournament_selection(population, fitness_scores))
# Perform crossover with crossover probability
offspring_list = []
for i in range(0, len(parents_list) - 1, 2):
# Apply crossover based on crossover probability
if random.random() < crossover_per:
offspring_1, offspring_2 = sbx_crossover(
parents_list[i], parents_list[i + 1], eta=1.5
)
else:
# If no crossover, children are exact copies of parents
offspring_1, offspring_2 = parents_list[i][:], parents_list[i + 1][:]
offspring_list.extend([offspring_1, offspring_2])
# Perform mutation
for i in range(len(offspring_list)):
if random.random() < mutation_per:
offspring_list[i] = gaussian_mutation(offspring_list[i])
# Combine elite individuals, parents and offspring
mixed_population = offspring_list + elite_individuals
# Filter solutions that meet the mission constraints
mixed_population = [
ind
for ind in mixed_population
if np.abs(ind[0] - date_departure_mjd) <= change_in_time_days
and np.abs(ind[1] - date_flyby_jupiter_mjd) <= change_in_time_days
and np.abs(ind[2] - date_flyby_saturn_mjd) <= change_in_time_days
and np.abs(ind[3] - date_flyby_uranus_mjd) <= change_in_time_days
and np.abs(ind[4] - date_flyby_neptune_mjd) <= change_in_time_days
]
# Ensure population is not empty (if all individuals were filtered out)
if not mixed_population:
mixed_population = (
offspring_list + elite_individuals
) # Restore full population if filtering removed all
# Ensure population size remains constant (selection 'n_population' best individuals)
fitness_scores = np.array([fitness(ind) for ind in mixed_population])
population = sort_population_by_fitness(
mixed_population, fitness_scores, n_population
)
# Store current population
population_per_generation.append(population)
end_time = time.time()
# print(end_time - start_time)
# Return the best-found solution
best_individual = population[0]
best_list = best_deltaV, best_fuel_mass, best_data = find_deltaV_mission(
*best_individual
)
print(
f"Generation {n_generations}: Best Delta-V = {best_deltaV:.6f} km/s, Fitness = {fitness(best_individual):.6f}"
)
print("\nFinal Optimised Solution:")
print(f"Departure Date (MJD): {best_individual[0]:.2f}")
print(f"Jupiter Flyby Date (MJD): {best_individual[1]:.2f}")
print(f"Saturn Flyby Date (MJD): {best_individual[2]:.2f}")
print(f"Uranus Flyby Date (MJD): {best_individual[3]:.2f}")
print(f"Neptune Flyby Date (MJD): {best_individual[4]:.2f}")
print(f"Total deltaV: {best_deltaV:.6f} km/s")
return best_individual, best_list, population_per_generation
In [130]:
# n_population = 2000
# n_elite = 5
# crossover_per = 0.9
# mutation_per = 0.35
# n_generations = 25
# n_population = 2500
# n_elite = 10
# crossover_per = 0.9
# mutation_per = 0.3
# n_generations = 300
n_population = 500
n_elite = 25
crossover_per = 0.9
mutation_per = 0.1
n_generations = 100
In [131]:
%%time
if global_flag:
# Start timing
start_time = time.time()
# Run the Genetic Algorithm
best_individual, best_mission_data, population_per_generation = (
run_genetic_algorithm(
n_population, n_elite, crossover_per, mutation_per, n_generations, []
)
)
print()
# End timing
end_time = time.time()
execution_time_GA = end_time - start_time
print(f"\nExecution Time: {execution_time_GA:.4f} seconds\n")
CPU times: user 1e+03 ns, sys: 0 ns, total: 1e+03 ns Wall time: 2.86 µs
In [132]:
if global_flag:
filename_GA = output_pickle.get_path("ga_results_V2.pkl")
results_GA = {
"Best Individual": best_individual,
"Best Mission Data": best_mission_data,
"Population Per Generation": population_per_generation,
"execution_time_GA": execution_time_GA,
}
with open(filename_GA, "wb") as f:
pickle.dump(results_GA, f)
In [133]:
filename_GA = output_pickle.get_path("ga_results_V2.pkl")
with open(filename_GA, "rb") as f:
data = pickle.load(f)
best_individual = data["Best Individual"]
best_mission_data = data["Best Mission Data"]
population_per_generation = data["Population Per Generation"]
execution_time_GA = data["execution_time_GA"]
if not global_flag:
for gen in range(len(population_per_generation)):
best_individual = population_per_generation[gen][0]
best_deltaV, best_fuel_mass, _ = find_deltaV_mission(*best_individual)
print(
f"Generation {gen}: Best Delta-V = {best_deltaV:.6f} km/s, Fitness = {fitness(best_individual):.6f}"
)
# Return the best-found solution
best_individual = population_per_generation[-1][0]
best_list = best_deltaV, best_fuel_mass, best_data = find_deltaV_mission(
*best_individual
)
print(
f"Generation {n_generations}: Best Delta-V = {best_deltaV:.6f} km/s, Fitness = {fitness(best_individual):.6f}"
)
print("\nFinal Optimised Solution:")
print(f"Departure Date (MJD): {best_individual[0]:.2f}")
print(f"Jupiter Flyby Date (MJD): {best_individual[1]:.2f}")
print(f"Saturn Flyby Date (MJD): {best_individual[2]:.2f}")
print(f"Total deltaV: {best_deltaV:.6f} km/s")
print(f"\nExecution Time: {execution_time_GA:.4f} seconds\n")
Generation 0: Best Delta-V = 14.928106 km/s, Fitness = 0.066988 Generation 1: Best Delta-V = 14.928106 km/s, Fitness = 0.066988 Generation 2: Best Delta-V = 12.171653 km/s, Fitness = 0.082158 Generation 3: Best Delta-V = 10.979255 km/s, Fitness = 0.091081 Generation 4: Best Delta-V = 10.830302 km/s, Fitness = 0.092334 Generation 5: Best Delta-V = 10.528690 km/s, Fitness = 0.094979 Generation 6: Best Delta-V = 10.400745 km/s, Fitness = 0.096147 Generation 7: Best Delta-V = 10.400745 km/s, Fitness = 0.096147 Generation 8: Best Delta-V = 10.254958 km/s, Fitness = 0.097514 Generation 9: Best Delta-V = 10.145931 km/s, Fitness = 0.098562 Generation 10: Best Delta-V = 9.978617 km/s, Fitness = 0.100214 Generation 11: Best Delta-V = 9.860274 km/s, Fitness = 0.101417 Generation 12: Best Delta-V = 9.763277 km/s, Fitness = 0.102425 Generation 13: Best Delta-V = 9.697365 km/s, Fitness = 0.103121 Generation 14: Best Delta-V = 9.694029 km/s, Fitness = 0.103156 Generation 15: Best Delta-V = 9.663997 km/s, Fitness = 0.103477 Generation 16: Best Delta-V = 9.607678 km/s, Fitness = 0.104083 Generation 17: Best Delta-V = 9.534269 km/s, Fitness = 0.104885 Generation 18: Best Delta-V = 9.488404 km/s, Fitness = 0.105392 Generation 19: Best Delta-V = 9.488404 km/s, Fitness = 0.105392 Generation 20: Best Delta-V = 9.425861 km/s, Fitness = 0.106091 Generation 21: Best Delta-V = 9.425861 km/s, Fitness = 0.106091 Generation 22: Best Delta-V = 9.419401 km/s, Fitness = 0.106164 Generation 23: Best Delta-V = 9.393977 km/s, Fitness = 0.106451 Generation 24: Best Delta-V = 9.390216 km/s, Fitness = 0.106494 Generation 25: Best Delta-V = 9.382785 km/s, Fitness = 0.106578 Generation 26: Best Delta-V = 9.380029 km/s, Fitness = 0.106609 Generation 27: Best Delta-V = 9.376201 km/s, Fitness = 0.106653 Generation 28: Best Delta-V = 9.363378 km/s, Fitness = 0.106799 Generation 29: Best Delta-V = 9.363378 km/s, Fitness = 0.106799 Generation 30: Best Delta-V = 9.362748 km/s, Fitness = 0.106806 Generation 31: Best Delta-V = 9.362017 km/s, Fitness = 0.106815 Generation 32: Best Delta-V = 9.362017 km/s, Fitness = 0.106815 Generation 33: Best Delta-V = 9.362017 km/s, Fitness = 0.106815 Generation 34: Best Delta-V = 9.361989 km/s, Fitness = 0.106815 Generation 35: Best Delta-V = 9.361989 km/s, Fitness = 0.106815 Generation 36: Best Delta-V = 9.361888 km/s, Fitness = 0.106816 Generation 37: Best Delta-V = 9.361798 km/s, Fitness = 0.106817 Generation 38: Best Delta-V = 9.361798 km/s, Fitness = 0.106817 Generation 39: Best Delta-V = 9.361688 km/s, Fitness = 0.106818 Generation 40: Best Delta-V = 9.361688 km/s, Fitness = 0.106818 Generation 41: Best Delta-V = 9.361688 km/s, Fitness = 0.106818 Generation 42: Best Delta-V = 9.361679 km/s, Fitness = 0.106818 Generation 43: Best Delta-V = 9.361640 km/s, Fitness = 0.106819 Generation 44: Best Delta-V = 9.361640 km/s, Fitness = 0.106819 Generation 45: Best Delta-V = 9.361634 km/s, Fitness = 0.106819 Generation 46: Best Delta-V = 9.361621 km/s, Fitness = 0.106819 Generation 47: Best Delta-V = 9.361613 km/s, Fitness = 0.106819 Generation 48: Best Delta-V = 9.361584 km/s, Fitness = 0.106820 Generation 49: Best Delta-V = 9.361548 km/s, Fitness = 0.106820 Generation 50: Best Delta-V = 9.361548 km/s, Fitness = 0.106820 Generation 51: Best Delta-V = 9.361548 km/s, Fitness = 0.106820 Generation 52: Best Delta-V = 9.361548 km/s, Fitness = 0.106820 Generation 53: Best Delta-V = 9.361548 km/s, Fitness = 0.106820 Generation 54: Best Delta-V = 9.361542 km/s, Fitness = 0.106820 Generation 55: Best Delta-V = 9.361535 km/s, Fitness = 0.106820 Generation 56: Best Delta-V = 9.361535 km/s, Fitness = 0.106820 Generation 57: Best Delta-V = 9.361527 km/s, Fitness = 0.106820 Generation 58: Best Delta-V = 9.361527 km/s, Fitness = 0.106820 Generation 59: Best Delta-V = 9.361523 km/s, Fitness = 0.106820 Generation 60: Best Delta-V = 9.361523 km/s, Fitness = 0.106820 Generation 61: Best Delta-V = 9.361523 km/s, Fitness = 0.106820 Generation 62: Best Delta-V = 9.361523 km/s, Fitness = 0.106820 Generation 63: Best Delta-V = 9.361523 km/s, Fitness = 0.106820 Generation 64: Best Delta-V = 9.361523 km/s, Fitness = 0.106820 Generation 65: Best Delta-V = 9.361522 km/s, Fitness = 0.106820 Generation 66: Best Delta-V = 9.361522 km/s, Fitness = 0.106820 Generation 67: Best Delta-V = 9.361522 km/s, Fitness = 0.106820 Generation 68: Best Delta-V = 9.361522 km/s, Fitness = 0.106820 Generation 69: Best Delta-V = 9.361522 km/s, Fitness = 0.106820 Generation 70: Best Delta-V = 9.361522 km/s, Fitness = 0.106820 Generation 71: Best Delta-V = 9.361522 km/s, Fitness = 0.106820 Generation 72: Best Delta-V = 9.361522 km/s, Fitness = 0.106820 Generation 73: Best Delta-V = 9.361522 km/s, Fitness = 0.106820 Generation 74: Best Delta-V = 9.361521 km/s, Fitness = 0.106820 Generation 75: Best Delta-V = 9.361521 km/s, Fitness = 0.106820 Generation 76: Best Delta-V = 9.361521 km/s, Fitness = 0.106820 Generation 77: Best Delta-V = 9.361521 km/s, Fitness = 0.106820 Generation 78: Best Delta-V = 9.361521 km/s, Fitness = 0.106820 Generation 79: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 80: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 81: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 82: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 83: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 84: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 85: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 86: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 87: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 88: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 89: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 90: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 91: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 92: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 93: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 94: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 95: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 96: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 97: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 98: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 99: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 100: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Generation 100: Best Delta-V = 9.361520 km/s, Fitness = 0.106820 Final Optimised Solution: Departure Date (MJD): 43388.81 Jupiter Flyby Date (MJD): 44116.31 Saturn Flyby Date (MJD): 44963.77 Total deltaV: 9.361520 km/s Execution Time: 634.4994 seconds
In [134]:
time.sleep(1) # Give the kernel a moment to catch up
In [135]:
best_solutions_per_generation = [element[0] for element in population_per_generation]
# best_solutions_per_generation, len(best_solutions_per_generation)
In [136]:
dates_REAL
Out[136]:
[43378.478599537164, 44063.93739583343, 44842.142326388974, 46666.74984953692, 47763.164305555634]
GA Result Analysis¶
In [138]:
deltaV_GA, fuel_mass_GA, data_dict_GA = find_deltaV_mission(*best_individual)
deltaV_GA, fuel_mass_GA, data_dict_GA
Out[138]:
(9.36152032789823,
816.1285632058713,
{'deltaV_injection (km/s)': 9.361520326659656,
'deltaV_jupiter (km/s)': 7.18996417958806e-10,
'rp_jupiter (km)': 2770004.045423409,
'vp_jupiter (km/s)': 11.901891880948948,
'vp_esc_jupiter (km/s)': 9.56500003203136,
'deflection_jupiter (deg)': 94.25132811528943,
'deltaV_saturn (km/s)': 2.1593926646801265e-10,
'rp_saturn (km)': 507478.59536726406,
'vp_saturn (km/s)': 15.441204304495137,
'vp_esc_saturn (km/s)': 12.22656521092557,
'deflection_saturn (deg)': 84.34962634943922,
'deltaV_uranus (km/s)': 3.03638003629203e-10,
'rp_uranus (km)': 168233.3340319893,
'vp_uranus (km/s)': 15.670785764613605,
'vp_esc_uranus (km/s)': 8.299385590522153,
'deflection_uranus (deg)': 20.627619471411386,
'tof (days)': 4716.798433988326})
In [139]:
date_departure_mjd_GA = best_individual[0]
date_flyby_jupiter_mjd_GA = best_individual[1]
date_flyby_saturn_mjd_GA = best_individual[2]
date_flyby_uranus_mjd_GA = best_individual[3]
date_flyby_neptune_mjd_GA = best_individual[4]
print(f"Difference for date_departure: {date_departure_mjd_GA - date_departure_mjd}")
print(
f"Difference for date_flyby_jupiter: {date_flyby_jupiter_mjd_GA - date_flyby_jupiter_mjd}"
)
print(
f"Difference for date_flyby_saturn: {date_flyby_saturn_mjd_GA - date_flyby_saturn_mjd}"
)
print(
f"Difference for date_flyby_uranus: {date_flyby_uranus_mjd_GA - date_flyby_uranus_mjd}"
)
print(
f"Difference for date_flyby_neptune: {date_flyby_neptune_mjd_GA - date_flyby_neptune_mjd}\n"
)
change_in_dates_list_GA = [
(date_departure_mjd_GA - date_departure_mjd),
(date_flyby_jupiter_mjd_GA - date_flyby_jupiter_mjd),
(date_flyby_saturn_mjd_GA - date_flyby_saturn_mjd),
(date_flyby_uranus_mjd_GA - date_flyby_uranus_mjd),
(date_flyby_neptune_mjd_GA - date_flyby_neptune_mjd),
]
change_in_dates_list_GA
Difference for date_departure: 10.330942864602548 Difference for date_flyby_jupiter: 52.3680094810843 Difference for date_flyby_saturn: 121.6322707088766 Difference for date_flyby_uranus: 23.88140782364644 Difference for date_flyby_neptune: 342.4436708344583
Out[139]:
[10.330942864602548, 52.3680094810843, 121.6322707088766, 23.88140782364644, 342.4436708344583]
In [140]:
print("Earth-Jupiter:", date_flyby_jupiter_mjd_GA - date_departure_mjd_GA, "days")
print("Jupiter-Saturn:", date_flyby_saturn_mjd_GA - date_flyby_jupiter_mjd_GA, "days")
print("Saturn-Uranus:", date_flyby_uranus_mjd_GA - date_flyby_saturn_mjd_GA, "days")
print("Uranus-Neptune:", date_flyby_neptune_mjd_GA - date_flyby_uranus_mjd_GA, "days")
Earth-Jupiter: 727.4958629127505 days Jupiter-Saturn: 847.4691917833334 days Saturn-Uranus: 1726.856660262718 days Uranus-Neptune: 1414.9767190295242 days
In [141]:
date_departure_GA = MJD_to_TT_calander_date(date_departure_mjd_GA)
date_flyby_jupiter_GA = MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA)
date_flyby_saturn_GA = MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA)
date_flyby_uranus_GA = MJD_to_TT_calander_date(date_flyby_uranus_mjd_GA)
date_flyby_neptune_GA = MJD_to_TT_calander_date(date_flyby_neptune_mjd_GA)
print("Date of Departure (Min):", MJD_to_TT_calander_date(date_departure_mjd_GA))
print(
"Date of Jupiter Flyby (Min):", MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA)
)
print("Date of Saturn Flyby (Min):", MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA))
print("Date of Uranus Flyby (Min):", MJD_to_TT_calander_date(date_flyby_uranus_mjd_GA))
print(
"Date of Neptune Flyby (Min):", MJD_to_TT_calander_date(date_flyby_neptune_mjd_GA)
)
[
date_departure_GA,
date_flyby_jupiter_GA,
date_flyby_saturn_GA,
date_flyby_uranus_GA,
date_flyby_neptune_GA,
]
Date of Departure (Min): 1977-09-02 19:25:44.464 Date of Jupiter Flyby (Min): 1979-08-31 07:19:47.019 Date of Saturn Flyby (Min): 1981-12-25 18:35:25.189 Date of Uranus Flyby (Min): 1986-09-17 15:09:00.636 Date of Neptune Flyby (Min): 1990-08-02 14:35:29.160
Out[141]:
['1977-09-02 19:25:44.464', '1979-08-31 07:19:47.019', '1981-12-25 18:35:25.189', '1986-09-17 15:09:00.636', '1990-08-02 14:35:29.160']
In [142]:
GA_list = [
date_departure_mjd_GA,
date_flyby_jupiter_mjd_GA,
date_flyby_saturn_mjd_GA,
date_flyby_uranus_mjd_GA,
date_flyby_neptune_mjd_GA,
]
GA_list
Out[142]:
[43388.80954240177, 44116.30540531452, 44963.77459709785, 46690.63125736057, 48105.60797639009]
Objective Function vs No. of Generations¶
In [144]:
generations = list(
range(len(best_solutions_per_generation))
) # X-axis (Generation numbers)
deltaV_values_GA = [] # Y-axis (Best Delta V per generation)
for solution in best_solutions_per_generation:
best_deltaV, best_fuel_mass, _ = find_deltaV_mission(
*solution
) # Compute Delta V for each best solution
deltaV_values_GA.append(best_deltaV)
# Plot the data
fig, ax = plt.subplots(figsize=(9, 6))
plt.plot(
generations[:],
deltaV_values_GA[:],
marker=".",
linestyle="--",
label="Best $\Delta V$ [km/s] per Generation",
)
plt.xlabel("No. of Generations")
plt.ylabel("Objective Function")
# plt.title(r"\textbf{Objective Function vs No. of Generations}", pad=10, fontsize=14)
plt.legend(loc="best")
# Improve grid and ticks
ax.grid(which="major", color="#DDDDDD", linewidth=0.8)
ax.grid(which="minor", color="#EEEEEE", linestyle=":", linewidth=0.5)
ax.minorticks_on()
# Customize tick parameters
ax.tick_params(axis="both", which="major", labelsize=12)
ax.tick_params(axis="both", which="minor", labelsize=10)
plt.tight_layout()
# Show plot
plt.show()
Fitness Function vs No. of Generations¶
In [146]:
fitness_values_deltaV = 1 / np.array(deltaV_values_GA)
fig, ax = plt.subplots(figsize=(9, 6))
plt.plot(
generations,
fitness_values_deltaV,
marker=".",
linestyle="--",
label="Fitness Function $(1/\Delta V)$ per Generation",
)
plt.xlabel("No. of Generations", fontsize=12)
plt.ylabel("Fitness Function", fontsize=12)
# plt.title(r"\textbf{Fitness Function vs No. of Generations}", pad=10, fontsize=14)
# Improve grid and ticks
ax.grid(which="major", color="#DDDDDD", linewidth=0.8)
ax.grid(which="minor", color="#EEEEEE", linestyle=":", linewidth=0.5)
ax.minorticks_on()
# Customise tick parameters
ax.tick_params(axis="both", which="major", labelsize=12)
ax.tick_params(axis="both", which="minor", labelsize=10)
plt.tight_layout()
output_figures.save_file(
"voyager_II_GA_fitness_function.png", dpi=500, bbox_inches="tight"
)
# Show plot
plt.show()
Figure saved to: Voyager 2 - Figures/voyager_II_GA_fitness_function.png
Objective & Fitness¶
In [148]:
fig, ax1 = plt.subplots(figsize=(9, 6))
# Plot objective value on primary axis (left)
color1 = "tab:red"
ax1.set_xlabel("No. of Generations", fontsize=12)
ax1.set_ylabel("Objective Function $(\\Delta V)$ [km/s]", fontsize=12, color=color1)
line1 = ax1.plot(
generations,
deltaV_values_GA,
marker=".",
linestyle="-",
linewidth=1,
color=color1,
label="Objective Function $(\\Delta V)$",
)
ax1.tick_params(axis="y", labelcolor=color1)
# Create secondary axis (right) for fitness value
ax2 = ax1.twinx()
color2 = "tab:blue"
ax2.set_ylabel("Fitness Function $(1/\\Delta V)$", fontsize=12, color=color2)
line2 = ax2.plot(
generations,
fitness_values_deltaV,
marker=".",
linestyle="-",
linewidth=1,
color=color2,
label="Fitness Function $(1/\\Delta V)$",
)
ax2.tick_params(axis="y", labelcolor=color2)
lines = line1 + line2
labels = [l.get_label() for l in lines]
ax1.legend(lines, labels, loc="upper right", bbox_to_anchor=(1, 0.94), fontsize=11)
ax1.grid(which="major", color="#DDDDDD", linewidth=0.8)
ax1.grid(which="minor", color="#EEEEEE", linestyle=":", linewidth=0.5)
ax1.minorticks_on()
# plt.title(r'\textbf{Fitness \& Objective Functions vs No. of Generations}',
# pad=10, fontsize=14)
# ax1.set_ylim(1/0.1068, 1/0.1028) # For objective function (ΔV)
# ax2.set_ylim(0.103, 0.1072) # For fitness function (1/ΔV)
ax1.set_xlim(-1, 51)
plt.tight_layout()
output_figures.save_file(
"voyager_II_GA_fitness_objective.png", dpi=500, bbox_inches="tight"
)
# Show plot
plt.show()
Figure saved to: Voyager 2 - Figures/voyager_II_GA_fitness_objective.png
In [149]:
gc.collect()
Out[149]:
15787
GA Results --> Fmin¶
In [151]:
%%time
if global_flag:
iterations = []
iteration_count = []
# Start timing
start_time = time.time()
# Run optimisation using `fmin` with manual bounds enforcement
optimal_dates = fmin(
objective_function,
GA_list,
disp=True,
xtol=1e-9,
ftol=1e-9,
maxiter=1e6,
maxfun=1e6,
)
iteration_count = len(iteration_count)
# Extract optimised dates
(
date_departure_mjd_GA_FMIN_COMBINED,
date_flyby_jupiter_mjd_GA_FMIN_COMBINED,
date_flyby_saturn_mjd_GA_FMIN_COMBINED,
date_flyby_uranus_mjd_GA_FMIN_COMBINED,
date_flyby_neptune_mjd_GA_FMIN_COMBINED,
) = optimal_dates
# Compute Delta-V and fuel mass at optimised dates
deltaV_GA_FMIN_COMBINED, fuel_mass_GA_FMIN_COMBINED, data_dict_GA_FMIN_COMBINED = (
find_deltaV_mission(*optimal_dates)
)
# Display results
print(f"\nMinimum Delta-V: {deltaV_GA_FMIN_COMBINED:.4f} km/s")
print(f"Optimal Departure Date (MJD): {date_departure_mjd_GA_FMIN_COMBINED:.2f}")
print(
f"Optimal Jupiter Flyby Date (MJD): {date_flyby_jupiter_mjd_GA_FMIN_COMBINED:.2f}"
)
print(
f"Optimal Saturn Flyby Date (MJD): {date_flyby_saturn_mjd_GA_FMIN_COMBINED:.2f}"
)
print(
f"Optimal Uranus Flyby Date (MJD): {date_flyby_uranus_mjd_GA_FMIN_COMBINED:.2f}"
)
print(
f"Optimal Neptune Flyby Date (MJD): {date_flyby_neptune_mjd_GA_FMIN_COMBINED:.2f}"
)
# Convert MJD dates to calendar dates
print("\nIn calendar dates:")
print(f"Departure: {MJD_to_TT_calander_date(date_departure_mjd_GA_FMIN_COMBINED)}")
print(
f"Jupiter flyby: {MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA_FMIN_COMBINED)}"
)
print(
f"Saturn flyby: {MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA_FMIN_COMBINED)}"
)
print(
f"Uranus flyby: {MJD_to_TT_calander_date(date_flyby_uranus_mjd_GA_FMIN_COMBINED)}"
)
print(
f"Neptune flyby: {MJD_to_TT_calander_date(date_flyby_neptune_mjd_GA_FMIN_COMBINED)}"
)
# End timing
execution_time_GA_FMIN_COMBINED = time.time() - start_time
print(f"\nExecution Time: {execution_time_GA_FMIN_COMBINED:.4f} seconds")
print(f"No. of Iterations: {iteration_count}")
CPU times: user 1e+03 ns, sys: 1 µs, total: 2 µs Wall time: 3.1 µs
In [152]:
if global_flag:
filename_brute = output_pickle.get_path("ga_+_fmin_results_V2.pkl")
results_brute = {
"optimal_dates": optimal_dates,
"execution_time_GA_FMIN_COMBINED": execution_time_GA_FMIN_COMBINED,
"iteration_count": iteration_count,
}
with open(filename_brute, "wb") as f:
pickle.dump(results_brute, f)
In [153]:
filename_brute = output_pickle.get_path("ga_+_fmin_results_V2.pkl")
with open(filename_brute, "rb") as f:
data = pickle.load(f)
optimal_dates = data["optimal_dates"]
execution_time_GA_FMIN_COMBINED = data["execution_time_GA_FMIN_COMBINED"]
iteration_count = data["iteration_count"]
if not global_flag:
# Extract optimised dates
(
date_departure_mjd_GA_FMIN_COMBINED,
date_flyby_jupiter_mjd_GA_FMIN_COMBINED,
date_flyby_saturn_mjd_GA_FMIN_COMBINED,
date_flyby_uranus_mjd_GA_FMIN_COMBINED,
date_flyby_neptune_mjd_GA_FMIN_COMBINED,
) = optimal_dates
# Compute Delta-V and fuel mass at optimised dates
deltaV_GA_FMIN_COMBINED, fuel_mass_GA_FMIN_COMBINED, data_dict_GA_FMIN_COMBINED = (
find_deltaV_mission(*optimal_dates)
)
# Display results
print(f"\nMinimum Delta-V: {deltaV_GA_FMIN_COMBINED:.4f} km/s")
print(f"Optimal Departure Date (MJD): {date_departure_mjd_GA_FMIN_COMBINED:.2f}")
print(
f"Optimal Jupiter Flyby Date (MJD): {date_flyby_jupiter_mjd_GA_FMIN_COMBINED:.2f}"
)
print(
f"Optimal Saturn Flyby Date (MJD): {date_flyby_saturn_mjd_GA_FMIN_COMBINED:.2f}"
)
print(
f"Optimal Uranus Flyby Date (MJD): {date_flyby_uranus_mjd_GA_FMIN_COMBINED:.2f}"
)
print(
f"Optimal Neptune Flyby Date (MJD): {date_flyby_neptune_mjd_GA_FMIN_COMBINED:.2f}"
)
# Convert MJD dates to calendar dates
print("\nIn calendar dates:")
print(f"Departure: {MJD_to_TT_calander_date(date_departure_mjd_GA_FMIN_COMBINED)}")
print(
f"Jupiter flyby: {MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA_FMIN_COMBINED)}"
)
print(
f"Saturn flyby: {MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA_FMIN_COMBINED)}"
)
print(
f"Uranus flyby: {MJD_to_TT_calander_date(date_flyby_uranus_mjd_GA_FMIN_COMBINED)}"
)
print(
f"Neptune flyby: {MJD_to_TT_calander_date(date_flyby_neptune_mjd_GA_FMIN_COMBINED)}"
)
print(f"\nExecution Time: {execution_time_GA_FMIN_COMBINED:.4f} seconds")
print(f"No. of Iterations: {iteration_count}")
Minimum Delta-V: 9.3580 km/s Optimal Departure Date (MJD): 43388.87 Optimal Jupiter Flyby Date (MJD): 44119.07 Optimal Saturn Flyby Date (MJD): 44970.95 Optimal Uranus Flyby Date (MJD): 46705.60 Optimal Neptune Flyby Date (MJD): 48127.84 In calendar dates: Departure: 1977-09-02 20:59:15.045 Jupiter flyby: 1979-09-03 01:38:02.079 Saturn flyby: 1982-01-01 22:47:12.250 Uranus flyby: 1986-10-02 14:31:04.214 Neptune flyby: 1990-08-24 20:04:15.847 Execution Time: 15.0133 seconds No. of Iterations: 3303
In [154]:
deltaV_GA_FMIN_COMBINED, fuel_mass_GA_FMIN_COMBINED, data_dict_GA_FMIN_COMBINED = (
find_deltaV_mission(*optimal_dates)
)
deltaV_GA_FMIN_COMBINED, fuel_mass_GA_FMIN_COMBINED, data_dict_GA_FMIN_COMBINED
Out[154]:
(9.358046430397152,
815.8954449993682,
{'deltaV_injection (km/s)': 9.358037063778704,
'deltaV_jupiter (km/s)': 8.953961387447862e-06,
'rp_jupiter (km)': 2800740.34130855,
'vp_jupiter (km/s)': 11.837743359432663,
'vp_esc_jupiter (km/s)': 9.51237042590936,
'deflection_jupiter (deg)': 94.15970492764589,
'deltaV_saturn (km/s)': 9.521272659185342e-13,
'rp_saturn (km)': 516021.4711421794,
'vp_saturn (km/s)': 15.314127824503991,
'vp_esc_saturn (km/s)': 12.124935773829252,
'deflection_saturn (deg)': 84.30292432886259,
'deltaV_uranus (km/s)': 4.126561083950264e-07,
'rp_uranus (km)': 171822.4661678232,
'vp_uranus (km/s)': 15.549357779710983,
'vp_esc_uranus (km/s)': 8.212246831191617,
'deflection_uranus (deg)': 20.479259426751486,
'tof (days)': 4738.961814830298})
In [155]:
execution_time_GA_FMIN_COMBINED_total = (
execution_time_GA + execution_time_GA_FMIN_COMBINED
)
execution_time_GA_FMIN_COMBINED_total
Out[155]:
649.5126950740814
In [156]:
print(
"Earth-Jupiter:",
date_flyby_jupiter_mjd_GA_FMIN_COMBINED - date_departure_mjd_GA_FMIN_COMBINED,
"days",
)
print(
"Jupiter-Saturn:",
date_flyby_saturn_mjd_GA_FMIN_COMBINED - date_flyby_jupiter_mjd_GA_FMIN_COMBINED,
"days",
)
print(
"Saturn-Uranus:",
date_flyby_uranus_mjd_GA_FMIN_COMBINED - date_flyby_saturn_mjd_GA_FMIN_COMBINED,
"days",
)
print(
"Uranus-Neptune:",
date_flyby_neptune_mjd_GA_FMIN_COMBINED - date_flyby_uranus_mjd_GA_FMIN_COMBINED,
"days",
)
Earth-Jupiter: 730.1935999267444 days Jupiter-Saturn: 851.8813677157596 days Saturn-Uranus: 1734.655462556082 days Uranus-Neptune: 1422.231384631712 days
In [157]:
change_in_dates_list_GA_FMIN_COMBINED = [
(date_departure_mjd_GA_FMIN_COMBINED - date_departure_mjd),
(date_flyby_jupiter_mjd_GA_FMIN_COMBINED - date_flyby_jupiter_mjd),
(date_flyby_saturn_mjd_GA_FMIN_COMBINED - date_flyby_saturn_mjd),
(date_flyby_uranus_mjd_GA_FMIN_COMBINED - date_flyby_uranus_mjd),
(date_flyby_neptune_mjd_GA_FMIN_COMBINED - date_flyby_neptune_mjd),
]
change_in_dates_list_GA_FMIN_COMBINED
Out[157]:
[10.395880153591861, 55.130683784067514, 128.807120944286, 38.85506035241997, 364.6719889654196]
In [158]:
date_departure_GA_FMIN_COMBINED = MJD_to_TT_calander_date(
date_departure_mjd_GA_FMIN_COMBINED
)
date_flyby_jupiter_GA_FMIN_COMBINED = MJD_to_TT_calander_date(
date_flyby_jupiter_mjd_GA_FMIN_COMBINED
)
date_flyby_saturn_GA_FMIN_COMBINED = MJD_to_TT_calander_date(
date_flyby_saturn_mjd_GA_FMIN_COMBINED
)
date_flyby_uranus_GA_FMIN_COMBINED = MJD_to_TT_calander_date(
date_flyby_uranus_mjd_GA_FMIN_COMBINED
)
date_flyby_neptune_GA_FMIN_COMBINED = MJD_to_TT_calander_date(
date_flyby_neptune_mjd_GA_FMIN_COMBINED
)
print(
"Date of Departure (Min):",
MJD_to_TT_calander_date(date_departure_mjd_GA_FMIN_COMBINED),
)
print(
"Date of Jupiter Flyby (Min):",
MJD_to_TT_calander_date(date_flyby_jupiter_mjd_GA_FMIN_COMBINED),
)
print(
"Date of Saturn Flyby (Min):",
MJD_to_TT_calander_date(date_flyby_saturn_mjd_GA_FMIN_COMBINED),
)
print(
"Date of Uranus Flyby (Min):",
MJD_to_TT_calander_date(date_flyby_uranus_mjd_GA_FMIN_COMBINED),
)
print(
"Date of Neptune Flyby (Min):",
MJD_to_TT_calander_date(date_flyby_neptune_mjd_GA_FMIN_COMBINED),
)
Date of Departure (Min): 1977-09-02 20:59:15.045 Date of Jupiter Flyby (Min): 1979-09-03 01:38:02.079 Date of Saturn Flyby (Min): 1982-01-01 22:47:12.250 Date of Uranus Flyby (Min): 1986-10-02 14:31:04.214 Date of Neptune Flyby (Min): 1990-08-24 20:04:15.847
In [159]:
[
date_departure_BRUTE_MIN,
date_flyby_jupiter_BRUTE_MIN,
date_flyby_saturn_BRUTE_MIN,
date_flyby_uranus_BRUTE_MIN,
date_flyby_neptune_BRUTE_MIN,
deltaV_BRUTE_MIN,
fuel_mass_BRUTE_MIN,
data_dict_BRUTE_MIN,
execution_time_BRUTE_FORCE,
]
Out[159]:
['1977-08-23 11:29:11.000',
'1979-07-09 22:29:51.000',
'1981-08-26 03:24:57.000',
'1986-01-28 04:16:55.571',
'1989-08-25 03:56:36.000',
9.965627065328311,
856.0135081245757,
{'deltaV_injection (km/s)': 9.777779726682427,
'deltaV_jupiter (km/s)': 0.03485226340972147,
'rp_jupiter (km)': 2220733.0033500125,
'vp_jupiter (km/s)': 13.264161047559059,
'vp_esc_jupiter (km/s)': 10.682600912576174,
'deflection_jupiter (deg)': 96.2960163290349,
'deltaV_saturn (km/s)': 0.0819711822100011,
'rp_saturn (km)': 385095.0255481553,
'vp_saturn (km/s)': 17.717014962912074,
'vp_esc_saturn (km/s)': 14.035548160975852,
'deflection_saturn (deg)': 85.03234833966007,
'deltaV_uranus (km/s)': 0.07102389302616174,
'rp_uranus (km)': 121875.94940717192,
'vp_uranus (km/s)': 17.652953215434174,
'vp_esc_uranus (km/s)': 9.750860929686262,
'deflection_uranus (deg)': 22.953523176610197,
'tof (days)': 4384.68570601847},
4812.23163485527]
Save Dates¶
In [161]:
if not testing_flag:
# Create a dictionary with all results
all_results = {
"real": [
date_departure,
date_flyby_jupiter,
date_flyby_saturn,
date_flyby_uranus,
date_flyby_neptune,
deltaV_REAL,
fuel_mass_REAL,
data_dict_REAL,
],
"brute_force": [
date_departure_BRUTE_MIN,
date_flyby_jupiter_BRUTE_MIN,
date_flyby_saturn_BRUTE_MIN,
date_flyby_uranus_BRUTE_MIN,
date_flyby_neptune_BRUTE_MIN,
deltaV_BRUTE_MIN,
fuel_mass_BRUTE_MIN,
data_dict_BRUTE_MIN,
change_in_dates_list_BRUTE_MIN,
execution_time_BRUTE_FORCE,
],
"fmin": [
date_departure_FMIN,
date_flyby_jupiter_FMIN,
date_flyby_saturn_FMIN,
date_flyby_uranus_FMIN,
date_flyby_neptune_FMIN,
deltaV_FMIN,
fuel_mass_FMIN,
data_dict_FMIN,
change_in_dates_list_FMIN,
execution_time_FMIN,
],
"ga": [
date_departure_GA,
date_flyby_jupiter_GA,
date_flyby_saturn_GA,
date_flyby_uranus_GA,
date_flyby_neptune_GA,
deltaV_GA,
fuel_mass_GA,
data_dict_GA,
change_in_dates_list_GA,
execution_time_GA,
],
"ga_fmin_combined": [
date_departure_GA_FMIN_COMBINED,
date_flyby_jupiter_GA_FMIN_COMBINED,
date_flyby_saturn_GA_FMIN_COMBINED,
date_flyby_uranus_GA_FMIN_COMBINED,
date_flyby_neptune_GA_FMIN_COMBINED,
deltaV_GA_FMIN_COMBINED,
fuel_mass_GA_FMIN_COMBINED,
data_dict_GA_FMIN_COMBINED,
change_in_dates_list_GA_FMIN_COMBINED,
execution_time_GA_FMIN_COMBINED,
],
}
# Save all results to a single pickle file
filename = output_pickle.get_path("voyager2_mission_dates.pkl")
with open(filename, "wb") as f:
pickle.dump(all_results, f)
print(f"Saved all mission dates to {filename}")
Saved all mission dates to Raw Results/voyager2_mission_dates.pkl
In [162]:
filename = output_pickle.get_path("voyager2_mission_dates.pkl")
# Read the data back from the pickle file
with open(filename, "rb") as f:
loaded_results = pickle.load(f)
print("\nLoaded data from pickle file:")
for method, data in loaded_results.items():
print(f"\n{method.upper()} mission dates:")
print(f" Departure: {data[0]}")
print(f" Jupiter flyby: {data[1]}")
print(f" Saturn flyby: {data[2]}")
print(f" Uranus flyby: {data[3]}")
print(f" Neptune flyby: {data[4]}")
print(f" Delta V: {data[5]:.6f} km/s")
print(f" Fuel Mass: {data[6]:.3f} kg")
# if len(data) > 8:
# print(f" Depature change in date: {data[-2][0]}")
# print(f" Jupiter change in date: {data[-2][1]}")
# print(f" Saturn change in date: {data[-2][2]}")
# print(f" Uranus change in date: {data[-2][3]}")
# print(f" Neptune change in date: {data[-2][4]}")
# print(f" Execution Time: {data[-1]:.3f} s")
print()
for i in data[7]:
print(f" {i}: {data[7][i]}")
print()
for method, data in loaded_results.items():
print(f"{method.upper()}:")
print(f" Delta V: {data[5]:.6f} km/s")
# Real Mission dates and data
(
date_departure,
date_flyby_jupiter,
date_flyby_saturn,
date_flyby_uranus,
date_flyby_neptune,
) = loaded_results["real"][:5]
deltaV_REAL, fuel_mass_REAL, data_dict_REAL = loaded_results["real"][5:]
# Brute Force Algorithm dates and data
(
date_departure_BRUTE_MIN,
date_flyby_jupiter_BRUTE_MIN,
date_flyby_saturn_BRUTE_MIN,
date_flyby_uranus_BRUTE_MIN,
date_flyby_neptune_BRUTE_MIN,
) = loaded_results["brute_force"][:5]
(
deltaV_BRUTE_MIN,
fuel_mass_BRUTE_MIN,
data_dict_BRUTE_MIN,
change_in_dates_list_BRUTE_MIN,
execution_time_BRUTE_FORCE,
) = loaded_results["brute_force"][5:]
# FMIN Algorithm dates
(
date_departure_FMIN,
date_flyby_jupiter_FMIN,
date_flyby_saturn_FMIN,
date_flyby_uranus_FMIN,
date_flyby_neptune_FMIN,
) = loaded_results["fmin"][:5]
(
deltaV_FMIN,
fuel_mass_FMIN,
data_dict_FMIN,
change_in_dates_list_FMIN,
execution_time_FMIN,
) = loaded_results["fmin"][5:]
# Genetic Algorithm dates
(
date_departure_GA,
date_flyby_jupiter_GA,
date_flyby_saturn_GA,
date_flyby_uranus_GA,
date_flyby_neptune_GA,
) = loaded_results["ga"][:5]
deltaV_GA, fuel_mass_GA, data_dict_GA, change_in_dates_list_GA, execution_time_GA = (
loaded_results["ga"][5:]
)
# Combined GA + FMIN dates
(
date_departure_GA_FMIN_COMBINED,
date_flyby_jupiter_GA_FMIN_COMBINED,
date_flyby_saturn_GA_FMIN_COMBINED,
date_flyby_uranus_GA_FMIN_COMBINED,
date_flyby_neptune_GA_FMIN_COMBINED,
) = loaded_results["ga_fmin_combined"][:5]
(
deltaV_GA_FMIN_COMBINED,
fuel_mass_GA_FMIN_COMBINED,
data_dict_GA_FMIN_COMBINED,
change_in_dates_list_GA_FMIN_COMBINED,
execution_time_GA_FMIN_COMBINED,
) = loaded_results["ga_fmin_combined"][5:]
execution_time_GA_FMIN_COMBINED_total = (
execution_time_GA + execution_time_GA_FMIN_COMBINED
)
print(
f"\n[GA + Fmin] Combined Execution Time: {execution_time_GA_FMIN_COMBINED_total:.3f} s"
)
Loaded data from pickle file: REAL mission dates: Departure: 1977-08-23 11:29:11 Jupiter flyby: 1979-07-09 22:29:51 Saturn flyby: 1981-08-26 03:24:57 Uranus flyby: 1986-08-24 17:59:47 Neptune flyby: 1989-08-25 03:56:36 Delta V: 14.928106 km/s Fuel Mass: 1139.000 kg deltaV_injection (km/s): 9.777779726682427 deltaV_jupiter (km/s): 0.03485226340972147 rp_jupiter (km): 2220733.0033500125 vp_jupiter (km/s): 13.264161047559059 vp_esc_jupiter (km/s): 10.682600912576174 deflection_jupiter (deg): 96.2960163290349 deltaV_saturn (km/s): 1.0288185212545642 rp_saturn (km): 403827.9428746656 vp_saturn (km/s): 17.457212514618856 vp_esc_saturn (km/s): 13.706139564574247 deflection_saturn (deg): 85.73086421478894 deltaV_uranus (km/s): 4.086655221591039 rp_uranus (km): 133605.1322173773 vp_uranus (km/s): 15.891238081241008 vp_esc_uranus (km/s): 9.313017043900006 deflection_uranus (deg): 22.164348557554597 tof (days): 4384.68570601847 BRUTE_FORCE mission dates: Departure: 1977-08-23 11:29:11.000 Jupiter flyby: 1979-07-09 22:29:51.000 Saturn flyby: 1981-08-26 03:24:57.000 Uranus flyby: 1986-01-28 04:16:55.571 Neptune flyby: 1989-08-25 03:56:36.000 Delta V: 9.965627 km/s Fuel Mass: 856.014 kg deltaV_injection (km/s): 9.777779726682427 deltaV_jupiter (km/s): 0.03485226340972147 rp_jupiter (km): 2220733.0033500125 vp_jupiter (km/s): 13.264161047559059 vp_esc_jupiter (km/s): 10.682600912576174 deflection_jupiter (deg): 96.2960163290349 deltaV_saturn (km/s): 0.0819711822100011 rp_saturn (km): 385095.0255481553 vp_saturn (km/s): 17.717014962912074 vp_esc_saturn (km/s): 14.035548160975852 deflection_saturn (deg): 85.03234833966007 deltaV_uranus (km/s): 0.07102389302616174 rp_uranus (km): 121875.94940717192 vp_uranus (km/s): 17.652953215434174 vp_esc_uranus (km/s): 9.750860929686262 deflection_uranus (deg): 22.953523176610197 tof (days): 4384.68570601847 FMIN mission dates: Departure: 1977-09-02 23:16:50.675 Jupiter flyby: 1979-09-03 05:06:38.181 Saturn flyby: 1982-01-02 11:28:33.046 Uranus flyby: 1986-10-03 18:21:11.480 Neptune flyby: 1990-08-25 03:56:35.928 Delta V: 9.371258 km/s Fuel Mass: 816.782 kg deltaV_injection (km/s): 9.35789550243903 deltaV_jupiter (km/s): 0.0009168183205279234 rp_jupiter (km): 2802517.371817503 vp_jupiter (km/s): 11.834271643601488 vp_esc_jupiter (km/s): 9.50935412829766 deflection_jupiter (deg): 94.14765367850407 deltaV_saturn (km/s): 2.6201263381153694e-11 rp_saturn (km): 516702.31495290494 vp_saturn (km/s): 15.304134784770053 vp_esc_saturn (km/s): 12.11694480065466 deflection_saturn (deg): 84.29926641857209 deltaV_uranus (km/s): 0.012445180557890012 rp_uranus (km): 171918.01181510068 vp_uranus (km/s): 15.542222697049406 vp_esc_uranus (km/s): 8.209964483111364 deflection_uranus (deg): 20.471397443524044 tof (days): 4739.194273761466 GA mission dates: Departure: 1977-09-02 19:25:44.464 Jupiter flyby: 1979-08-31 07:19:47.019 Saturn flyby: 1981-12-25 18:35:25.189 Uranus flyby: 1986-09-17 15:09:00.636 Neptune flyby: 1990-08-02 14:35:29.160 Delta V: 9.361520 km/s Fuel Mass: 816.129 kg deltaV_injection (km/s): 9.361520326659656 deltaV_jupiter (km/s): 7.18996417958806e-10 rp_jupiter (km): 2770004.045423409 vp_jupiter (km/s): 11.901891880948948 vp_esc_jupiter (km/s): 9.56500003203136 deflection_jupiter (deg): 94.25132811528943 deltaV_saturn (km/s): 2.1593926646801265e-10 rp_saturn (km): 507478.59536726406 vp_saturn (km/s): 15.441204304495137 vp_esc_saturn (km/s): 12.22656521092557 deflection_saturn (deg): 84.34962634943922 deltaV_uranus (km/s): 3.03638003629203e-10 rp_uranus (km): 168233.3340319893 vp_uranus (km/s): 15.670785764613605 vp_esc_uranus (km/s): 8.299385590522153 deflection_uranus (deg): 20.627619471411386 tof (days): 4716.798433988326 GA_FMIN_COMBINED mission dates: Departure: 1977-09-02 20:59:15.045 Jupiter flyby: 1979-09-03 01:38:02.079 Saturn flyby: 1982-01-01 22:47:12.250 Uranus flyby: 1986-10-02 14:31:04.214 Neptune flyby: 1990-08-24 20:04:15.847 Delta V: 9.358046 km/s Fuel Mass: 815.895 kg deltaV_injection (km/s): 9.358037063778704 deltaV_jupiter (km/s): 8.953961387447862e-06 rp_jupiter (km): 2800740.34130855 vp_jupiter (km/s): 11.837743359432663 vp_esc_jupiter (km/s): 9.51237042590936 deflection_jupiter (deg): 94.15970492764589 deltaV_saturn (km/s): 9.521272659185342e-13 rp_saturn (km): 516021.4711421794 vp_saturn (km/s): 15.314127824503991 vp_esc_saturn (km/s): 12.124935773829252 deflection_saturn (deg): 84.30292432886259 deltaV_uranus (km/s): 4.126561083950264e-07 rp_uranus (km): 171822.4661678232 vp_uranus (km/s): 15.549357779710983 vp_esc_uranus (km/s): 8.212246831191617 deflection_uranus (deg): 20.479259426751486 tof (days): 4738.961814830298 REAL: Delta V: 14.928106 km/s BRUTE_FORCE: Delta V: 9.965627 km/s FMIN: Delta V: 9.371258 km/s GA: Delta V: 9.361520 km/s GA_FMIN_COMBINED: Delta V: 9.358046 km/s [GA + Fmin] Combined Execution Time: 649.513 s
Graph¶
In [164]:
def plot_VoyagerII(
date_departure,
date_flyby_jupiter,
date_flyby_saturn,
date_flyby_uranus,
date_flyby_neptune,
colour,
label_suffix,
):
# Convert input dates to Astropy Time in TDB scale
date_departure = Time(date_departure, scale="tt").tdb
date_flyby_jupiter = Time(date_flyby_jupiter, scale="tt").tdb
date_flyby_saturn = Time(date_flyby_saturn, scale="tt").tdb
date_flyby_uranus = Time(date_flyby_uranus, scale="tt").tdb
date_flyby_neptune = Time(date_flyby_neptune, scale="tt").tdb
date_arrival = date_flyby_neptune
# Define planetary ephemerides
earth = Ephem.from_body(
Earth, time_range(date_departure, end=date_arrival, periods=500)
) # , plane=Planes.EARTH_ECLIPTIC)
ss_earth = Orbit.from_ephem(Sun, earth, date_departure)
jupiter = Ephem.from_body(
Jupiter, time_range(date_departure, end=date_arrival, periods=500)
) # , plane=Planes.EARTH_ECLIPTIC)
ss_jupiter = Orbit.from_ephem(Sun, jupiter, date_flyby_jupiter)
saturn = Ephem.from_body(
Saturn, time_range(date_departure, end=date_arrival, periods=500)
) # , plane=Planes.EARTH_ECLIPTIC)
ss_saturn = Orbit.from_ephem(Sun, saturn, date_flyby_saturn)
uranus = Ephem.from_body(
Uranus, time_range(date_departure, end=date_arrival, periods=500)
) # , plane=Planes.EARTH_ECLIPTIC)
ss_uranus = Orbit.from_ephem(Sun, uranus, date_flyby_uranus)
neptune = Ephem.from_body(
Neptune, time_range(date_departure, end=date_arrival, periods=500)
) # , plane=Planes.EARTH_ECLIPTIC)
ss_neptune = Orbit.from_ephem(Sun, neptune, date_arrival)
# Solving for maneuver to Jupiter
man_flyby_jupiter = Maneuver.lambert(ss_earth, ss_jupiter)
ic1 = ss_earth.apply_maneuver(man_flyby_jupiter)
ic1_end = ic1.propagate(date_flyby_jupiter)
# Solving for maneuver to Saturn
man_flyby_saturn = Maneuver.lambert(ic1_end, ss_saturn)
ic2 = ic1_end.apply_maneuver(man_flyby_saturn)
ic2_end = ic2.propagate(date_flyby_saturn)
# Solving for maneuver to Uranus
man_uranus = Maneuver.lambert(ic2_end, ss_uranus)
ic3 = ic2_end.apply_maneuver(man_uranus)
ic3_end = ic3.propagate(date_flyby_uranus)
# Solving for maneuver to Neptune
man_neptune = Maneuver.lambert(ic3_end, ss_neptune)
ic4 = ic3_end.apply_maneuver(man_neptune)
ic4_end = ic4.propagate(date_arrival)
# Plot Earth's position at departure
plotter.plot_body_orbit(
Earth, date_departure, label=f"Earth Departure {label_suffix}", trail=True
)
# Plot cruise to Jupiter
plotter.plot_maneuver(
ss_earth,
man_flyby_jupiter,
label=f"Cruise to Jupiter {label_suffix}",
color=colour,
)
# Plot Jupiter's position at flyby
plotter.plot_body_orbit(
Jupiter, date_flyby_jupiter, label=f"Jupiter Flyby {label_suffix}", trail=True
)
# Plot cruise to Saturn
plotter.plot_maneuver(
ic1_end,
man_flyby_saturn,
label=f"Cruise to Saturn {label_suffix}",
color=colour,
)
# Plot Saturn's position at arrival
plotter.plot_body_orbit(
Saturn, date_flyby_saturn, label=f"Saturn Arrival {label_suffix}", trail=True
)
# Plot cruise to Uranus
plotter.plot_maneuver(
ic2_end, man_uranus, label=f"Cruise to Uranus {label_suffix}", color=colour
)
# Plot Uranus's position at arrival
plotter.plot_body_orbit(
Uranus, date_flyby_uranus, label=f"Uranus Arrival {label_suffix}", trail=True
)
# Plot cruise to Neptune
plotter.plot_maneuver(
ic3_end, man_neptune, label=f"Cruise to Neptune {label_suffix}", color=colour
)
# Plot Neptune's position at arrival
plotter.plot_body_orbit(
Neptune, date_flyby_neptune, label=f"Neptune Arrival {label_suffix}", trail=True
)
2-Axis¶
In [166]:
set_font_sizes(14)
# Define different line styles for each plot
colours = {
"Real Mission": "red",
"Brute Force Search": "cyan",
"Fmin": "magenta",
"Genetic Algorithm": "gold",
"Hybrid GA": "lime",
}
# Distinct line styles to prevent visual merging
line_styles = {
"Real Mission": "-",
"Brute Force Search": "-",
"Fmin": "--",
"Genetic Algorithm": "--",
"Hybrid GA": "--",
}
# Define different line widths for clarity
line_widths = {
"Real Mission": 2,
"Brute Force Search": 2,
"Fmin": 3.5,
"Genetic Algorithm": 3,
"Hybrid GA": 1.75,
}
plt.rc("text", usetex=True)
plt.rc("font", family="serif")
# Create figure and axis
fig = plt.figure(figsize=(12, 10))
ax = fig.add_subplot(111)
plotter = StaticOrbitPlotter(ax)
# Plot both missions on the same graph
plot_VoyagerII(
date_departure,
date_flyby_jupiter,
date_flyby_saturn,
date_flyby_uranus,
date_flyby_neptune,
colours["Real Mission"],
label_suffix="Real Mission",
)
plot_VoyagerII(
date_departure_BRUTE_MIN,
date_flyby_jupiter_BRUTE_MIN,
date_flyby_saturn_BRUTE_MIN,
date_flyby_uranus_BRUTE_MIN,
date_flyby_neptune_BRUTE_MIN,
colours["Brute Force Search"],
label_suffix="Brute Force Search",
)
plot_VoyagerII(
date_departure_FMIN,
date_flyby_jupiter_FMIN,
date_flyby_saturn_FMIN,
date_flyby_uranus_FMIN,
date_flyby_neptune_FMIN,
colours["Fmin"],
label_suffix="Fmin",
)
plot_VoyagerII(
date_departure_GA,
date_flyby_jupiter_GA,
date_flyby_saturn_GA,
date_flyby_uranus_GA,
date_flyby_neptune_GA,
colours["Genetic Algorithm"],
label_suffix="Genetic Algorithm",
)
plot_VoyagerII(
date_departure_GA_FMIN_COMBINED,
date_flyby_jupiter_GA_FMIN_COMBINED,
date_flyby_saturn_GA_FMIN_COMBINED,
date_flyby_uranus_GA_FMIN_COMBINED,
date_flyby_neptune_GA_FMIN_COMBINED,
colours["Hybrid GA"],
label_suffix="Hybrid GA",
)
# Set fixed limits for primary axis
limit = 5.5e9
ax.set_xlim(-limit, limit)
ax.set_ylim(-limit, limit)
ax.set_aspect("equal") # Square aspect ratio
# Conversion factor (1 AU = 149,597,870.691 km)
km_to_AU = 1 / 149597870.691
# Create second set of axes for AU scale
# These won't be linked to the primary axes in the problematic way
ax_top = fig.add_axes(ax.get_position(), frameon=False)
ax_top.xaxis.tick_top()
ax_top.yaxis.tick_right()
ax_top.set_xlim(ax.get_xlim()[0] * km_to_AU, ax.get_xlim()[1] * km_to_AU)
ax_top.set_ylim(ax.get_ylim()[0] * km_to_AU, ax.get_ylim()[1] * km_to_AU)
# Improve axis labels with LaTeX formatting
ax.set_xlabel(
r"$x$ (km)",
) # fontsize=12)
ax.set_ylabel(
r"$y$ (km)",
) # fontsize=12)
ax_top.set_xlabel(
r"$x$ (AU)",
) # fontsize=12)
ax_top.set_ylabel(
r"$y$ (AU)",
) # fontsize=12)
ax_top.xaxis.set_label_position("top")
ax_top.yaxis.set_label_position("right")
# # Add grid
ax.grid(which="major", color="#CCCCCC", linewidth=0.5, linestyle="-", alpha=0.8)
ax.grid(which="minor", color="#DDDDDD", linewidth=0.25, linestyle=":", alpha=0.5)
ax.minorticks_on()
ax_top.grid(False)
ax_top.minorticks_on()
# Customise tick parameters
ax.tick_params(axis="both", which="major", labelsize=10, width=1, length=5, pad=5)
ax.tick_params(axis="both", which="minor", labelsize=8, width=0.5, length=3)
ax_top.tick_params(
axis="both", which="major", labelsize=10, width=1, length=5, pad=5, colors="#444444"
)
ax_top.tick_params(
axis="both", which="minor", labelsize=8, width=0.5, length=3, colors="#444444"
)
# Improved title with LaTeX formatting
# plt.title(r"\textbf{Voyager II Mission Trajectory Optimisation}", pad=10)
# Function to update AU axes when primary axes change
def update_au_axes(event=None):
ax_top.set_xlim(ax.get_xlim()[0] * km_to_AU, ax.get_xlim()[1] * km_to_AU)
ax_top.set_ylim(ax.get_ylim()[0] * km_to_AU, ax.get_ylim()[1] * km_to_AU)
fig.canvas.draw_idle()
# Connect the update function to the 'xlim_changed' and 'ylim_changed' events
ax.callbacks.connect("xlim_changed", update_au_axes)
ax.callbacks.connect("ylim_changed", update_au_axes)
# Apply line styles
for line in ax.get_lines():
label = line.get_label()
for key in line_styles:
if key in label:
line.set_linestyle(line_styles[key])
line.set_linewidth(line_widths[key])
# Create custom legend
handles = [
Line2D(
[0],
[0],
color=colours["Real Mission"],
ls=line_styles["Real Mission"],
lw=line_widths["Real Mission"],
label="Real Mission",
),
Line2D(
[0],
[0],
color=colours["Brute Force Search"],
ls=line_styles["Brute Force Search"],
lw=line_widths["Brute Force Search"],
label="Brute Force Search",
),
Line2D(
[0],
[0],
color=colours["Fmin"],
ls=line_styles["Fmin"],
lw=line_widths["Fmin"],
label="Fmin", # r"\texttt{fmin}",
),
Line2D(
[0],
[0],
color=colours["Genetic Algorithm"],
ls=line_styles["Genetic Algorithm"],
lw=line_widths["Genetic Algorithm"],
label="Genetic Algorithm (GA)",
),
Line2D(
[0],
[0],
color=colours["Hybrid GA"],
ls=line_styles["Hybrid GA"],
lw=line_widths["Hybrid GA"],
label=r"Hybrid GA (GA + Fmin)",
),
]
# Sun
ax.plot(
0,
0,
"o",
markersize=8,
markerfacecolor="#FFFF00",
markeredgecolor="#FF8C00",
markeredgewidth=1.5,
zorder=100,
)
# Add legend with improved styling
legend = ax.legend(
handles=handles,
loc="upper left",
ncol=1,
frameon=True,
edgecolor="black",
# fontsize=10, # title=r"\textbf{Optimisation Methods}",
# title_fontsize=11,
)
legend.get_frame().set_alpha(0.9)
legend.get_frame().set_facecolor("#F9F9F9")
# Planet colors
planet_colors = {
"Earth": "#3366FF", # Blue
"Jupiter": "#FF9933", # Orange
"Saturn": "#FFCC33", # Yellow
}
# Use the new label positions provided
planet_label_positions = {
"Earth": (-0.35e9, -0.5e9),
"Jupiter": (0.9e9, 0.5e9),
"Saturn": (-0.3e9, 1.75e9),
"Uranus": (-2.95e9, 1.4e9),
"Neptune": (-4.925e9, 0.5e9),
}
# Add each planet marker and label
for planet, position in planet_label_positions.items():
x, y = position
# Use the absolute positions for labels instead of offsets
label_x, label_y = planet_label_positions[planet]
# Improved annotation style
bbox_props = dict(boxstyle="round,pad=0.4", fc="white", ec="black", alpha=0.85)
ax.annotate(
planet,
xy=(x, y), # Position of the planet
xytext=(label_x, label_y), # Absolute position for the label
textcoords="data", # Changed from 'offset points' to 'data'
bbox=bbox_props,
#
fontsize=11,
fontweight="bold",
zorder=101, # Ensure labels appear on top of everything
) # Add arrow connecting label to planet
output_figures.save_file(
"voyager_II_trajectory_km_AU.png", dpi=500, bbox_inches="tight"
)
# Show plot
plt.show()
plt.close()
Figure saved to: Voyager 2 - Figures/voyager_II_trajectory_km_AU.png
2-Axis Zoomed¶
In [168]:
# Define different line widths for clarity
# line_widths = {
# "Real Mission": 2.5,
# "Brute Force Search": 1.5,
# "Fmin": 3.5,
# "Genetic Algorithm": 3,
# "Hybrid GA": 1.75,
# }
set_font_sizes(13.5)
fig, ax = plt.subplots(figsize=(9, 7))
plotter = StaticOrbitPlotter(ax)
# Plot both missions on the same graph
plot_VoyagerII(
date_departure,
date_flyby_jupiter,
date_flyby_saturn,
date_flyby_uranus,
date_flyby_neptune,
colours["Real Mission"],
label_suffix="Real Mission",
)
plot_VoyagerII(
date_departure_BRUTE_MIN,
date_flyby_jupiter_BRUTE_MIN,
date_flyby_saturn_BRUTE_MIN,
date_flyby_uranus_BRUTE_MIN,
date_flyby_neptune_BRUTE_MIN,
colours["Brute Force Search"],
label_suffix="Brute Force Search",
)
plot_VoyagerII(
date_departure_FMIN,
date_flyby_jupiter_FMIN,
date_flyby_saturn_FMIN,
date_flyby_uranus_FMIN,
date_flyby_neptune_FMIN,
colours["Fmin"],
label_suffix="Fmin",
)
plot_VoyagerII(
date_departure_GA,
date_flyby_jupiter_GA,
date_flyby_saturn_GA,
date_flyby_uranus_GA,
date_flyby_neptune_GA,
colours["Genetic Algorithm"],
label_suffix="Genetic Algorithm",
)
plot_VoyagerII(
date_departure_GA_FMIN_COMBINED,
date_flyby_jupiter_GA_FMIN_COMBINED,
date_flyby_saturn_GA_FMIN_COMBINED,
date_flyby_uranus_GA_FMIN_COMBINED,
date_flyby_neptune_GA_FMIN_COMBINED,
colours["Hybrid GA"],
label_suffix="Hybrid GA",
)
# Set fixed limits for primary axis
plt.xlim(-5.1e9, 2e9)
plt.ylim(-1.75e9, 2.5e9)
ax.set_aspect("equal") # Square aspect ratio
# Conversion factor (1 AU = 149,597,870.691 km)
km_to_AU = 1 / 149597870.691
# Create second set of axes for AU scale
# These won't be linked to the primary axes in the problematic way
ax_top = fig.add_axes(ax.get_position(), frameon=False)
ax_top.xaxis.tick_top()
ax_top.yaxis.tick_right()
ax_top.set_xlim(ax.get_xlim()[0] * km_to_AU, ax.get_xlim()[1] * km_to_AU)
ax_top.set_ylim(ax.get_ylim()[0] * km_to_AU, ax.get_ylim()[1] * km_to_AU)
# Improve axis labels with LaTeX formatting
ax.set_xlabel(
r"$x$ (km)",
) # fontsize=12)
ax.set_ylabel(
r"$y$ (km)",
) # fontsize=12)
ax_top.set_xlabel(
r"$x$ (AU)",
) # fontsize=12)
ax_top.set_ylabel(
r"$y$ (AU)",
) # fontsize=12)
ax_top.xaxis.set_label_position("top")
ax_top.yaxis.set_label_position("right")
# # Add grid
ax.grid(which="major", color="#CCCCCC", linewidth=0.5, linestyle="-", alpha=0.8)
ax.grid(which="minor", color="#DDDDDD", linewidth=0.25, linestyle=":", alpha=0.5)
ax.minorticks_on()
ax_top.grid(False)
ax_top.minorticks_on()
# Customise tick parameters
ax.tick_params(axis="both", which="major", labelsize=10, width=1, length=5, pad=5)
ax.tick_params(axis="both", which="minor", labelsize=8, width=0.5, length=3)
ax_top.tick_params(
axis="both", which="major", labelsize=10, width=1, length=5, pad=5, colors="#444444"
)
ax_top.tick_params(
axis="both", which="minor", labelsize=8, width=0.5, length=3, colors="#444444"
)
# Improved title with LaTeX formatting
# plt.title(r"\textbf{Voyager II Mission Trajectory Optimisation}", pad=10)
# Function to update AU axes when primary axes change
def update_au_axes(event=None):
ax_top.set_xlim(ax.get_xlim()[0] * km_to_AU, ax.get_xlim()[1] * km_to_AU)
ax_top.set_ylim(ax.get_ylim()[0] * km_to_AU, ax.get_ylim()[1] * km_to_AU)
fig.canvas.draw_idle()
# Connect the update function to the 'xlim_changed' and 'ylim_changed' events
ax.callbacks.connect("xlim_changed", update_au_axes)
ax.callbacks.connect("ylim_changed", update_au_axes)
# Apply line styles
for line in ax.get_lines():
label = line.get_label()
for key in line_styles:
if key in label:
line.set_linestyle(line_styles[key])
line.set_linewidth(line_widths[key])
# Create custom legend
handles = [
Line2D(
[0],
[0],
color=colours["Real Mission"],
ls=line_styles["Real Mission"],
lw=line_widths["Real Mission"],
label="Real Mission",
),
Line2D(
[0],
[0],
color=colours["Brute Force Search"],
ls=line_styles["Brute Force Search"],
lw=line_widths["Brute Force Search"],
label="Brute Force Search",
),
Line2D(
[0],
[0],
color=colours["Fmin"],
ls=line_styles["Fmin"],
lw=line_widths["Fmin"],
label="Fmin", # r"\texttt{fmin}",
),
Line2D(
[0],
[0],
color=colours["Genetic Algorithm"],
ls=line_styles["Genetic Algorithm"],
lw=line_widths["Genetic Algorithm"],
label="Genetic Algorithm (GA)",
),
Line2D(
[0],
[0],
color=colours["Hybrid GA"],
ls=line_styles["Hybrid GA"],
lw=line_widths["Hybrid GA"],
label=r"Hybrid GA (GA + Fmin)",
),
]
# Sun
ax.plot(
0,
0,
"o",
markersize=8,
markerfacecolor="#FFFF00",
markeredgecolor="#FF8C00",
markeredgewidth=1.5,
zorder=100,
)
# Add legend with improved styling
legend = ax.legend(
handles=handles,
loc="upper left",
ncol=1,
frameon=True,
edgecolor="black",
# fontsize=10, # title=r"\textbf{Optimisation Methods}",
# title_fontsize=11,
)
legend.get_frame().set_alpha(0.9)
legend.get_frame().set_facecolor("#F9F9F9")
# Planet colors
planet_colors = {
"Earth": "#3366FF", # Blue
"Jupiter": "#FF9933", # Orange
"Saturn": "#FFCC33", # Yellow
}
# Use the new label positions provided
planet_label_positions = {
"Earth": (-0.35e9, -0.5e9),
"Jupiter": (0.9e9, 0.5e9),
"Saturn": (-0.3e9, 1.75e9),
"Uranus": (-2.7e9, 1.4e9),
"Neptune": (-5e9, 0.4e9),
}
# Add each planet marker and label
for planet, position in planet_label_positions.items():
x, y = position
# Use the absolute positions for labels instead of offsets
label_x, label_y = planet_label_positions[planet]
# Improved annotation style
bbox_props = dict(boxstyle="round,pad=0.4", fc="white", ec="black", alpha=0.85)
ax.annotate(
planet,
xy=(x, y), # Position of the planet
xytext=(label_x, label_y), # Absolute position for the label
textcoords="data", # Changed from 'offset points' to 'data'
bbox=bbox_props,
# fontsize=11,
fontweight="bold",
zorder=101, # Ensure labels appear on top of everything
) # Add arrow connecting label to planet
output_figures.save_file("voyager_II_trajectory_zoom.png", dpi=500, bbox_inches="tight")
# Show plot
plt.show()
plt.close()
Figure saved to: Voyager 2 - Figures/voyager_II_trajectory_zoom.png
Genetic Algorithm (GA) Animation¶
Orbital Trajectory per Generation¶
In [171]:
if not testing_flag:
# Ensure real mission is first in the best solutions list
best_solutions_per_generation, len(best_solutions_per_generation)
In [172]:
if not testing_flag:
fig, ax = plt.subplots(figsize=(9, 7))
plt.close(fig) # This prevents the initial empty figure from showing
plotter = StaticOrbitPlotter(ax)
# Store previously plotted trajectories
plotted_trajectories = []
# Select a colour map (viridis, plasma, coolwarm, rainbow, etc.)
cmap = plt.cm.rainbow # Change this for different color gradients
# Flag to track if we've plotted the real mission trajectory
real_mission_plotted = False
# Function to update the plot at each generation
def update(frame):
global real_mission_plotted
# Plot the real mission trajectory only once at the first frame
if frame == 0:
# Plot real mission in red
real_traj = plot_VoyagerII(
date_departure,
date_flyby_jupiter,
date_flyby_saturn,
date_flyby_uranus,
date_flyby_neptune,
"red",
label_suffix="(Real Mission)",
)
plotted_trajectories.append(real_traj)
real_mission_plotted = True
else:
# Get best trajectory for this generation
best_traj = best_solutions_per_generation[frame - 1]
date_departure_gen, date_flyby_jupiter_gen, date_flyby_saturn_gen, date_flyby_uranus_gen, date_flyby_neptune_gen = (
best_traj
)
# Convert MJD to calendar date
date_departure_gen = MJD_to_TT_calander_date(date_departure_gen)
date_flyby_jupiter_gen = MJD_to_TT_calander_date(date_flyby_jupiter_gen)
date_flyby_saturn_gen = MJD_to_TT_calander_date(date_flyby_saturn_gen)
date_flyby_uranus_gen = MJD_to_TT_calander_date(date_flyby_uranus_gen)
date_flyby_neptune_gen = MJD_to_TT_calander_date(date_flyby_neptune_gen)
# Normalise color based on generation number
norm_color = cmap(frame / len(best_solutions_per_generation))
# Plot new trajectory with gradient color
traj_plot = plot_VoyagerII(
date_departure_gen,
date_flyby_jupiter_gen,
date_flyby_saturn_gen,
date_flyby_uranus_gen,
date_flyby_neptune_gen,
norm_color,
label_suffix=f"(Gen {frame-1})",
)
# Store the trajectory plot for reference
plotted_trajectories.append(traj_plot)
# Set fixed limits for primary axis
ax.set_xlim(-5.1e9, 2e9)
ax.set_ylim(-1.75e9, 2.5e9)
ax.set_aspect("equal") # Square aspect ratio
# Hide the legend during animation to prevent clutter
ax.legend().set_visible(False)
if frame == 0:
ax.set_title(r"\textbf{Voyager II GA Optimisation - Real Mission}")
else:
# Set title dynamically
ax.set_title(
rf"\textbf{{Voyager II GA Optimisation - Generation {frame-1}}}"
)
return plotted_trajectories
# Improve grid and ticks
ax.grid(which="major", color="#DDDDDD", linewidth=0.8)
ax.grid(which="minor", color="#EEEEEE", linestyle=":", linewidth=0.5)
ax.minorticks_on()
# Customise tick parameters
ax.tick_params(axis="both", which="major", )#labelsize=12)
ax.tick_params(axis="both", which="minor", )#labelsize=10)
plt.tight_layout()
# Create the animation
ani = animation.FuncAnimation(
fig,
update,
frames=len(best_solutions_per_generation) + 1,
repeat=False,
blit=False,
)
# Save the animation as a GIF or MP4
# ani.save("voyager_I_GA_optimisation.mp4", writer="ffmpeg", fps=5)#, dpi=200, bitrate=-1)
# ani.save("voyager_I_GA_optimisation.gif", writer="pillow", fps=10, dpi=200) # Save as GIF
<Figure size 640x480 with 0 Axes>
In [173]:
if not testing_flag:
# Convert animation to HTML
html_video = ani.to_jshtml()
# Display animation inline
display(HTML(html_video))